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We give a factorization formula for the e+e^ thrust distribution dcr/dr with r = 1 — T based 
on soft-colhnear effective theory. The result is applicable for all r, i.e. in the peak, tail, and far- 
tail regions. The formula includes 0{al) fixed-order QCD results, resummation of singular partonic 
ai In*" (r) /r terms with N'^LL accuracy, hadronization effects from fitting a universal nonperturbative 
soft function defined in field theory, bottom quark mass effects, QED corrections, and the dominant 
top mass dependent terms from the axial anomaly. We do not rely on Monte Carlo generators 
to determine nonperturbative effects since they are not compatible with higher order perturbative 
analyses. Instead our treatment is based on fitting nonperturbative matrix elements in field theory, 
which are moments Qi of a nonperturbative soft function. We present a global analysis of all available 
thrust data measured at center-of-mass energies Q = 35 to 207 GeV in the tail region, where a two 
parameter fit to as{mz) and the first moment Qi suffices. We use a short distance scheme to 
define Sli, called the R-gap scheme, thus ensuring that the perturbative dcr/dr does not suffer 
from an C'(Aqcd) renormalon ambiguity. We find as{mz) = 0.1135 ± (0.0002)oxpt ± (0.0005)hadr ± 
(0.0009)pcrt, with x^/dof = 0.91, where the displayed 1-sigma errors are the total experimental 
error, the hadronization uncertainty, and the perturbative theory uncertainty, respectively. The 
hadronization uncertainty in as is significantly decreased compared to earlier analyses by our two 
parameter fit, which determines fii — 0.323 GeV with 16% uncertainty. 



I. INTRODUCTION 

A traditional method for testing the theory of strong 
interactions (QCD) at high-precision is the analysis of 
jet cross sections at e~ colliders. Event shape distri- 
butions play a special role as they have been extensively 
measured with small experimental uncertainties at LEP 
and earlier e+ colliders, and are theoretically clean 
and accessible to high-order perturbative computations. 
They have been frequently used to make precise determi- 
nations of the strong coupling a^, see e.g. Ref. [l| for a 
review. One of the most frequently studied event shape 
variables is thrust Q], 



T 



(1) 



where the sum i is over all final-state hadrons with mo- 
menta Pi. The unit vector t that maximizes the right- 
hand side (RHS) of Eq. ([T]) defines the thrust axis. We 
will use the more convenient variable r = 1 — T. For 
the production of a pair of massless quarks at tree level 
dcr/dr oc ^(t), so the measured distribution for r > 
involves gluon radiation and is sensitive to the value of 
ag. The thrust value of an event measures how much it 
resembles two jets. For r values close to zero the event 
has two narrow, pencil-like, back-to-back jets, carrying 
about half the center-of-mass (cm.) energy into each of 
the two hemispheres defined by the plane orthogonal to 
t. For r close to the multijet endpoint 1/2, the event has 
an isotropic multi-particle final state containing a large 
number of low-energy jets. 



On the theoretical side, for r < 1/3 the dynamics 
is governed by three different scales. The hard scale 
fiH — Q is set by the e~^e~ cm. energy Q. The jet 
scale, jij ~ Q\/t is the typical momentum transverse to 
t of the particles within each of the two hemispheres, or 
the jet invariant mass scale if all energetic particles in a 
hemisphere are grouped into a jet. There is also uniform 
soft radiation with energy ns — Qt, called the soft scale. 
The physical description of the thrust distribution can 
be divided into three regions. 



peak region: 
tail region: 
far-tail region: 



r - 2Aqcd/Q , 
2AQCD/Q«r <l/3, 
1/3 < r < 1/2. 



(2) 



In the peak region the hard. jet. and soft scales are 
Q- vQAqcd, and Aqcd, and the distribution shows a 
strongly peaked maximum. Theoretically, since r ^ 1 
one needs to sum large (double) logarithms, {a^ ln'^r)/r, 
and account for the fact that ^.s — Aqcd, so dcr/dr is 
affected at leading order by a nonperturbative distribu- 
tion. We call this distribution the nonperturbative soft 
function. The tail region is populated predominantly by 
broader dijets and 3-jet events. Here the three scales 
are still well separated and one still needs to sum loga- 
rithms, but now /i5 ^ Aqcd, so soft radiation can be 
described by perturbation theory and a series of power 
correction parameters fi^. Finally, the far-tail region is 
populated by multijet events. Here the distinction of 
the three scales becomes meaningless, and accurate pre- 
dictions can be made with fixed-order perturbation the- 
ory supplemented with power corrections. The transition 
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sum logs 


power corrections 


as{mz) 


Ref. 


[221 


no 


Monte Carlo (MC) 


0.1240 ±0.0034* 


Ref. 


[201 


N^LL 


uncertainty from MC 0.1172 ± 0.0021* 


Ref. 


[23j 


NLL 


effective coupling 


0.1164 ± 0.0028* 






model 




Ref. 


m 


NLL 


Monte Carlo 


0.1172 ± 0.0051** 


Ref. 


[25j 


NLL 


Monte Carlo 


0.1224 ± 0.0039* 



TABLE I: Recent thrust analyses which use the 0{al) fixed- 
order results. The theoretical component of the errors were 
determined as indicated, by either: * the error band method, 
** variation of the renormalization scale ^, or * by a simul- 
taneous fit to as{mz) and ao (see text for more details). 
The analyses of Refs. [20l . [23| used thrust data only, while 
Refs. [53, 113, [2^ employed six different event shapes. 



to this region must be handled carefully since including 
a summation of (a-jln''r)/r terms in this region spoils 
the cancellations that take place at fixed order multijct 
thresholds, and hence would induce uncertainties that are 
significantly larger than those of the fixed-order results. 

Recently two very important achievements were made 
improving the theoretical description of event shape dis- 
tributions in e^e~ annihilation. First, in the work of 
Gehrmann et al. in Refs. and Weinzierl in Refs. 
the full set of 0{a'^) contributions to the 2-, 3- and 4-jet 
final states were determined. These results were made 
available in the program package EERAD3 0. Second, 
soft-collinear effective theory (SCET) [7l-[TTj provides a 
systematic framework to treat nonperturbative correc- 
tions and to factorize and compute hard, collinear 
and soft contributions for jet production to all orders in 
buil ding on earlier all orders QCD factoriza- 
tion results [iMil- The SCET framework allows for 
the summation of large logarithms at higher orders, as 
demonstrated by the analytic calculation of the thrust 
distribution at N^LL order by Becher and Schwartz in 
Ref. [lO].^ In contrast, the classic exponentiation tech- 
niques of Ref. [m for event shapes have so far only been 
carried out to NLL order. Also, the anomalous dimen- 
sions in SCET relevant for thrust are valid over pertur- 
bative momentum scales, and there are no Landau pole 
ambiguities in the resummation at any order. In addi- 
tion, as we will discuss in the body of our paper, SCET 
provides a rigorous framework for including perturbative 
and nonperturbative contributions, which can be used 
to connect power corrections in factorization theorems 
to those in an operator expansion for thrust moments. 
Moreover it provides a simple method to simultaneously 
treat the peak, tail, and far-tail regions. 

Several determinations of in the tail region have 
been carried out incorporating the fixed-order 0{a^) re- 



^ The calculation of Ref. [20l | also revealed a numerical problem at 
small T in the initial fixed-order results of Refs. 0,0]. 



suits, which we have collected in Tab. HI They differ on 
which event shape data has been used for the fits, on 
the accuracy of the partonic resummation of logarithms 
in the theory formula, the approach for nonperturbative 
hadronization effects, and how the theory errors are esti- 
mated. It is instructive to compare the analyses by Dis- 
scrtori ct al. [2^ and by Becher and Schwartz , 
which both used the error band method [26| to determine 
theoretical uncertainties. The improved convergence and 
reduced theoretical uncertainty for as(rnz) obtained by 
Becher and Schwartz indicates that the summation of 
logarithms beyond NLL order level is important. Both 
the analyses by Dissertori et al. and Becher and Schwartz 
are limited by the fact that they used Monte Carlo (MC) 
generators to estimate the size of nonperturbative cor- 
rections. 

The use of e"'"e~ MC generators to estimate power cor- 
rections is problematic since the partonic contributions 
arc based on LL parton showers with at most one-loop 
matrix elements, complemented by hadronization models 
below the shower cutoff that are not derived from QCD. 
The parameters of these models have been tuned to LEP 
data, and thus unavoidably encode both nonperturba- 
tive effects as well as higher order perturbative correc- 
tions. Hence, one must worry about double counting, and 
this makes MC generators unreliable for estimating non- 
perturbative corrections in higher order LEP analyses. 
Moreover, purely perturbative results for event shapes in 
the MS scheme such as those in Refs. [1, [13, [12, [Ig suf- 
fer from infrared effects known as infrared renormalons 
(see Ref. [l^l for a review of the early literature). These 
infrared effects arise because fluctuations from large an- 
gle soft radiation down to arbitrarily small momenta are 
included in the MS perturbative series and can cause 
unphysical large corrections already in low-order per- 
turbative QCD results. On the other hand, the hard 
shower cutoff protects the parton level MC from infrared 
renormalons. Hence one cannot rigorously combine MC 
hadronization effects with strict perturbative MS results. 
From the two points raised above, we conclude that the 
as('7iz) results obtained in Refs. [lO, [lH, [1^ contain a 
systematic theoretical error from nonperturbative effects 
that can be quite sizeable. We emphasize that this crit- 
icism also applies in part to the numerous earlier event 
shape analyses which estimated nonperturbative correc- 
tions using MC generators, see Ref. [l[ for a review. 

The presence of Aqcd/(Qt) power corrections in 
dcr/dr have been discussed in earlier literature [28l - [33 |. 
where it has been argued that the leading effect is a shift 
in the thrust distribution, r — > t — 2A/Q with A ~ Aqcd- 
The analyses with the 0{al) results, that discuss nonper- 
turbative effects in the thrust tail region without relying 
on MC generators are: Ref. [2^ which examined a. 1/Q 
power correction in a simple soft function model, but due 
to the large induced uncertainty on as (rri z) d oes not use 
it for their final error analysis, and Ref. [23 which uses 
the effective coupling model. 

For the most accurate data at Q = mz the change to 
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FIG. 1: Plot of h'(T)/h{T), the slope of ln[(l/r)dcr/dr], com- 
puted from experimental data at Q = mz- The derivative 
is computed using the central difference with neighboring ex- 
perimental bins. 



the extracted as{mz) from including the leading power 
correction can be quite significant, at the 10% level. Wc 
can derive this estimate by a simple calculation. We first 
write the cross section with a shift due to the power cor- 
rection, {l/a)da/dT = /i(r — 2A/Q). Assuming h is pro- 
portional to as, and expanding for A/Q <^ 1. one can 
easily derive that the change in the value extracted 
from data due to the existence of the power correction is 



6as _ 2A /i'(t) 
as ~ Q h{T) 



(3) 



The expression in Eq. ^ gives a scaling estimate for 
the fractional change in as from an analysis with the 
power correction compared to one without, using data 
at T. To the extent that the assumptions stated above 
are realized, the slope factor h'{T)/h{T) should be con- 
stant. In Fig. [T]wc show the slope factor computed from 
experimental data ai Q ~ mz- The figure shows that 
the slope factor is indeed only weakly depending on r in 
the tail region and we get h'{T)/h{T) ~ —14 ± 4. The 
remaining visible variation in r is related to subleading 
nonperturbativc and higher power as effects that are not 
accounted for in our simple scaling estimate. For a QCD 
power correction of natural size, A — 0.3 GcV, Eq. (|3]) 
gives 5as/as — —(9 ± 3)% for Q = mz- The magnitude 
of this effect makes it important to treat power correc- 
tions as accurately as possible in a fit to thrust data. We 
will show in later sections of this work that the relative 
downward shift in the fitted as(mz) due to nonperturba- 
tivc effects is indeed at the level of the scaling estimate 
of -10%. 

In the NLL/C'(a^) analysis by Davison and Web- 
ber [2^ the nonperturbativc effects are incorporated 
through a power correction a^ which is fit together with 
as to the experimental data. The power correction is 
formulated from the low-scale effective coupling model of 
Ref. [35|, which modifies as{fJ.) below fj, = = 2 GcV, 



coupling model correctly predicts the Q dependence of 
the lea ding nonperturbativc power correction in factor- 
ization |l8l. l35j . This model also induces a subtraction of 
perturbative contributions below the momentum cutoff 
fii (based on the running coupling approximation) and 
thus removes infrared rcnormalon ambiguities.^ However 
the model is not based on factorization, and hence this 
treatment of nonperturbativc corrections is not system- 
atically improvable. It is therefore not easy to judge the 
corresponding uncertainty. Another problem of the effec- 
tive coupling model is that its subtractions involve large 
logs, \n{^i/Q), which are not resummed. This effects the 
Q dependence in the interplay between perturbative and 
nonperturbativc effects. 

In this paper we extend the event shape formalism to 
resolve the theoretical difficulties mentioned above. Our 
results are formulated in the SCET framework, and hence 
are rigorous predictions of QCD. The formula wc derive 
has a N'^LL order summation of logarithms for the par- 
tonic singular a| In*^ (t)/t terms, and 0{a^) fixed-order 
contributions for the partonic nonsingular terms. Our 
theoretical improvements beyond earlier work include: 

• A factorization formula that can be simultaneously 
applied to data in the peak and the tail regions of 
the thrust distribution and for multiple cm. ener- 
gies Q, as well as being consistent with the multijet 
thresholds in the far-tail region. 

• In the factorization formula a nonperturbativc 
soft function defined from field theory is imple- 
mented using the method of Ref. [3 to incorporate 
hadronization effects. To achieve independence of 
a particular analytic ansatz in the peak region, the 
nonperturbativc part of the soft function uses a lin- 
ear combination of orthogonal basis functions that 
converge quickly for confined functions [ssj . 

• In the tail region the leading power correction to 
da/dr is determined by a nonperturbativc param- 
eter f2i that appears through a factorization the- 
orem for the singular distribution, fii is a field 
theory matrix element of an operator, and is also 
related to the first moment of the nonperturbativc 
soft function. In the tail region the effects of J7i 
hadronization corrections are included for the non- 
singular corrections that arc kincmatically sublead- 
ing in the dijct limit, based on theoretical consis- 
tency with the far-tail region. 

• Defining the matrix element Hi in MS, the pertur- 
bative cross section suffer from an O(Aqcd) renor- 



^ Another thrust analysis where infrared renormalon contributions 
have been removed from the partonic contributions is by Gardi 



and defines ao as the average value of the coupling be- 
tween ^ = and /i/. It is important that the effective 



and Rathsman in Refs. 1361 . [33, which used a principal value 
prescription for the inverse Borel transformation of the thrust 
distribution. Their analysis was prior to the new O(of) fixed- 
order computations, and hence was not included in Tabic [l] 
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malon. In our analysis this renormalon is removed 
by using an R-gap scheme for the definition of 
rii [iBl- This scheme choice induces subtractions 
on the leading power MS cross section which si- 
multaneously remove the renormalon there. Large 
logarithms in the subtractions are summed to all 
orders in as using R-evolution equations given in 
Refs. 

• Finite bottom quark mass corrections are included 
using a factorization theorem for event shap es in- 
volving massive quarks, derived in Refs. [3|4l|. 

• QED corrections at NNLL order arc incorporated, 
counting acm ~ ce^- This includes QED Sudakov 
effects, final state radiation, and QED/QCD renor- 
malization group interference. 

• The 3-loop finite term of the quark form factor 
in MS is extracted using the results of Ref. (i^ . 
and is included in our analysis. 

• The most important corrections from the axial 
anomaly are included. The anomaly modifies 
the axial- vector current contributions at 0{a1) by 
terms involving the top quark mass. 

Electroweak effects from virtual W and Z loops mostly 
effect the normalization of the cross section and so their 
dominant contribution drops out of {l/a)da/dT (43ll4^. 
These corrections arc not included in our analysis. 

For the numerical analyses carried out in this work 
we have created within our collaboration two completely 
independent codes. One code within Mathematica [45| 
implements the theoretical expressions exactly as given 
in this paper, and one code is based on theoretical formu- 
lae in Fourier space and realized as a fast Fortran code 
suitable for parallelized runs on computer clusters. These 
two codes agree for the thrust distribution at the level of 
10-6. 

While the resulting theoretical code can be used for all 
values of t, in this paper we focus our numerical analysis 
on a global fit of e+e" thrust data in the tail region, for 
cm. energies Q between 35 and 207 GeV, to determine 
cks("^z)-'^ Our global fit exhibits consistency across all 
available data sets, and reduces the overall experimental 
uncertainty. For a single Q we find a strong correlation 
between the effect of as{mz) and i7i on the cross sec- 
tion. This degeneracy is broken by fitting data at mul- 
tiple Qs. The hadronization uncertainty is significantly 
decreased by our simultaneous global fit to as (mz) and 
f2i. To estimate the perturbative uncertainty in the fit we 
use a random scan in a 12-dimensional theory parameter 
space. This space includes 6 parameters for /i-variation, 
3 parameters for theory uncertainties related to the finite 



^ Throughout this paper wc use the MS scheme for as with five 
hght flavors. 



statistics of the numerical fixed-order results, one param- 
eter for the unknown 4-loop cusp anomalous dimension, 
and two parameters for unknown constants in the per- 
turbative 3-loop jet and 3-loop soft functions. The scan 
yields a more conservative theory error than the error 
band method [l^. Despite this we are able to achieve 
smaller perturbative uncertainties than earlier analyses 
due to our removal of the C'(Aqcd) renormalon and the 
inclusion of . We also analyze in detail the dependence 
of the fit results on the range in r used in the fit. 

The outline of the paper is as follows. In Sec. |ll] we 
introduce the theoretical framework and discuss the var- 
ious theoretical ingredients in our final dcr/dr formula. 
In Sec. mil we present the profile functions which allow 
us to simultaneously treat multiple t regions, and dis- 
cuss the 6 parameters used for //-variation in the analy- 
sis of the perturbative uncertainty. In Sec. IIVI we review 
the parametrization of the nonperturbative function. In 
Sec. |V] we discuss the normalization of our distributions 
and compare results at different orders in perturbation 
theory for: fixed-order results, adding the log resumma- 
tion, adding the nonperturbative corrections, and adding 
the renormalon subtractions. In Sec. IVII we discuss the 
experimental data and the fit procedure. Our results 
for as{rnz) and the soft function moment fii from the 
global fit are presented in Sec. IVIIi including a discus- 
sion of the theory errors. In Sec. IVIIII we use our tail fit 
results to make predictions in the far-tail and peak re- 
gions, and compare with data. Cross checks on our code 
are discussed in Sec. IIXI including using it to re pro duce 
the earlier lower precision fits of Dissertori et al. [22| and 
Becher and Schwartz [23] . Section |X] contains our con- 
clusions and outlook, including prospects for future im- 
provements based on the universality of the parameter 
r^i. The analytic theoretical expressions that went into 
our analysis for massless quarks and QCD effects are pre- 
sented in condensed form in Appendix 1X1 In Appendix [B] 
we use the operator product expansion for the soft func- 
tion in the tail region, discussing uniqueness and deriving 
an all order relation for the Wilson coefficient of fii. In 
Appendix [C] we use an OPE for the first moment of the 
thrust distribution to show that it involves the same fii 
at lowest order. Readers most interested in our numerical 
results can skip directly to Sections |Vl] and IVIII 



II. FORMALISM 

A. Overview 

The factorization formula we use for the fits to the 
experimental thrust data is 

l\i ( do's , do-ns , Ado-(,\ / k\ 

d7=y^^ (,d7+^+^j(,"-Qj 

X Sr''{k-2^{R,^is)) +0(a^as^^) . (4) 
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Here dag/dr contains the singular partonic QCD correc- 
tions ai [In {t)/t]+ and al S{t) with the standard plus- 
functions as defined in Eq. (jA17p . It also contains the sin- 
gular partonic QED corrections depending on acm which 
are discussed in Sec. Ill HI This dds/dr term accounts for 
matrix element corrections and the resummation of In t 
terms within the SCET formalism up to N^LL order, 
which we discuss in Sec. IH CI Our definition of N'^LL, 
N'^LL', and other orders is discussed in detail in Sec. lHBl 
(see also Tab.|ll|. 

The term dtTns/dr, which we call the nonsingular 
partonic distribution, contains the thrust distribution 
in strict fixed-order expansion with the singular terms 
oc In'^ (r)/T subtracted to avoid double counting. The 
most singular terms in dans/dr scale as In*" r for r — > 0.** 
Our implementation of nonsingular terms is discussed in 
detail in Sec. IHEI 

Finally, Adah/dr contains corrections to the singular 
and nonsingular cross sections due to the finite mass 
of the bottom quark. The 6-mass corrections are im- 
plemented as a difference of the massive and massless 
cross sections computed at NNLL order as discussed in 
Scc.lHGl 

The function 5™°'^ that is convoluted with these par- 
tonic cross sections in Eq. ^ describes the nonpertur- 
bativc effects from soft gluons including large angle soft 
radiation [3 EB] • The definition of S^^°'^ also depends on 
the hemisphere prescription inherent to the thrust vari- 
able. This is a hadronic function that enters in a uni- 
versal way for both massless and massive cross sections, 
and is independent of the value of Q. The universality 
of 5;."°^ in Eq. (g]) follows from the leading power thrust 
factorization theorem [13, [H, [13: and the thrust fac- 
torization theorem for massive quarks in Refs. l4l| . 
Our treatment of the convolution of S™°'^ with diTns/dr 
yields a consistent treatment of multijet thresholds and 
the leading power correction to the operator expansion 
for the first moment of thrust. Details of our implemen- 
tation of power corrections and nonperturbative correc- 
tions are discussed in Sec. IH D] and Sec lIVI The function 
^mod normalized to unity and can be determined from 
experimental data. Its form depends on a gap parame- 
ter A and additional moment parameters fli which arc 
discussed below. 

The factorization formula given in Eq. ([4]) can be ap- 
plied simultaneously in the peak, tail, and the far-tail 
regions of Eq. ([2]), i.e. for all r values. In the peak re- 
gion d(7ns/dr is significantly smaller than dCTg/dr, and 
the full analytic form of the soft nonperturbative func- 
tion S™°'^{k) is relevant to determine the r-distribution 
since — Aqcd- Because /ijj ^ IJi-j ^ fJ^s, the sum- 
mation of logarithms of r is also crucial to achieve an 



accurate theoretical description. 

For much of the tail region the summation of In r terms 
remains important, although this is no longer the case 
when we reach t ~ 1/3. Likewise, the dominance of the 
singular partonic contributions remains as long as t < 
1/3, but the nonsingular terms become more important 
for increasing r (see Fig. [7] below). Near r ~ 1/3 the 
nonsingular terms become equal in size to the singular 
terms with opposite sign. Since 3> Aqcd in the tail 
region the effects of S*™""^ can be parameterized in terms 
of the moments 

= / (I) - 2A) , (5) 

where VIq ~ 1 since S*™""^ is normalized. Their impor- 
tance is determined by Vti/^Qry as discussed in Sec. Ill D[ 
so the first moment fii parameterizes the dominant 
power correction and higher moments provide increas- 
ingly smaller corrections. The first moment is defined 
by 

f^i ^ A + ^ (0 1 tr (0)r„ (0) id y„t (0)y ! (O) | O) , (6) 

where F,J(0) = Pexp(i(7 dsn ■ A{ns)), is similar 
but in the 3 representation, and we trace over color. Here 

id = 9{in ■ d—in ■ d) in ■ d + 9{in ■ d~in ■ d)in ■ d, (7) 

is a derivative operator^ involving light-like vectors n = 
(l,t) and n — (1,— t). fli is the field theory analog 
of the parameter ao employed in the low-scale effective 
coupling approach to power corrections. Since the renor- 
malon subtractions depend on a cutoff scale R and the 
rcnormalization scale /xs, all moments r2i(i?,/xs) as well 
as A{R, fis) are scale and scheme dependent quantities. 
The scheme we use to define f2i(i?, /ig) is described in 
Sec. Ill Fl In our fit to experimental data we use the R- 
gap scheme, and extract the first moment at a reference 
scale Ra = /xa = 2 GeV, i.e. we use A(i?A,MA) and 
hence Qi = rii(i?A,MA)- In the factorization theorem 
the gap appears evaluated at A(i?, /xg) and the scales 
(i?, /is) are connected to the reference scales (i?A,MA) 
using rcnormalization group equations. 

Finally, in the far-tail region t ~ 0.3 the singu- 
lar and the nonsingular partonic contributions dag/dr 
and d(Tns/dT become nearly equal with opposite signs, 
exhibiting a strong cancellation. This is due to the 
strong suppression of the fixed-order distribution in the 
three- and four-jet endpoint regions at r > 1/3 in fixed- 
order perturbation theory. In this region the summa- 
tion of logarithms of r must be switched off to avoid 



* For d(Tns/dT the resummation of In r terms is currently unknown. 
These terms could be determined with subleading factorization 
theorems in SCET. 



^ Note that id is defined in the cm. frame of the colliding e+e". 
One may also write id = fdr] e~^^^£Tiri) where £t(v) measures 
the sum of absolute transverse momenta at a given rapidity rj 
with respect to the thrust axis i llSl |47|| . 



6 



messing up this cancellation. Here our Eq. (|4]) re- 
duces to the pure fixed-order partonic thrust distribu- 
tion supplemented with power corrections coming from 
the convolution with the soft function. All three re- 
gions are smoothly joined together in Eq. Q. The 
proper summation (or non-summation) of logarithms 
is achieved through r-dependcnt renormalization scales, 
/.tj(r), fisir), and R{t) which wc call profile functions. 
They arc discussed in detail in Sec. IIIIl 

In the following subsections various ingredients of the 
factorization formula of Eq. ^ are presented in more de- 
tail. Compact results for the corresponding analytic ex- 
pressions for masslcss quarks in QCD are given in App.EI 
In Sees. Ill Gj and HTTll we describe how finite bottom mass 
and QED corrections are included in our analysis. The 
full formulae for these corrections will be presented in a 
future publication. 



B. Order Counting 

In the classic order counting used for fits to event 
shape distributions it is common to separately quote or- 
ders for the summation of logarithms and the fixed-order 
matching contributions. For fixed-order contributions 
the 0{as) contributions are called LO, the 0{al) contri- 
butions are called NLO, etc. This counting is motivated 
from the fact that at tree level the fixed-order thrust 
distribution vanishes for r > 0. For the summation one 
refers to LL (leading-log) summation if the one-loop cusp 
anomalous dimension is used to sum the double Sudakov 
logs, and NLL (next-to-leading-log) if the two-loop cusp 
and the one-loop non-cusp anomalous dimension terms 
are also included. 

In our analysis the summation orders (LL, NLL, ...) 
match the classical language. For the fixed-order contri- 
butions we account for the tree level 6{t) in LL and NLL, 
and we include 0{as) corrections at NLL' and NNLL, 
etc, as shown in Tab. Hlk . In SCET the summation can 
be carried out at both NNLL and N^LL The corre- 
sponding loop orders for the anomalous dimensions are 
also shown in Tab. UBi. Within SCET the summation of 
logarithms is achieved by renormalization group evolu- 
tion and the fixed-order corrections enter as scries eval- 
uated at each of the transition scales fin, fJ-j, and fis 
which we refer to as matching or matrix element correc- 
tions. The logs in the singular thrust cross section expo- 
nentiate to all orders if we use y, the Fourier-transformed 
variable to r. The orders we consider correspond to sum- 
ming the terms 
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oo 
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k=0 



N^LL 



clear the structure of the large logs that arc summed at 
each order. 

The nonsingular counting in Tab. ITik for the fixed- 
order series in dans/dr must be the same as for the 
matching and matrix element corrections to ensure that 
we exactly reproduce the fixed-order cross section when 
the rcsummcd result is expanded. Since the relative im- 
portance of the log resummation and the nonsingular 
terms varies depending on the r-region, we also consider 
an alternative "primed" counting scheme. In the primed 
counting all series for fixed-order quantities are included 
to one higher order in as- In this counting scheme the 
0{al) fixed-order results occur in N'^LL', which is the 
order we use for our final analysis. 

Also shown in Tab. Hlk arc columns for the fixed-order 
gap subtractions 6 — 6{R,^), and the gap anomalous 
dimensions 7^' . These terms are required to remove 
the leading ©(Aqcd) renormalon from the perturbative 
corrections, while still maintaining the same level of log 
resummation for terms in the cross section. The resum- 
mation of these large logarithms is missing in the recent 
analysis of Ref. [231 and is discussed further in Sec. Ill Fl 

A crucial aspect of our analysis is the inclusion of power 
corrections in a rigorous manner through field theoretic 
techniques. In the effective theory there are several types 
of power corrections which arise from the possible ratios 
of the scales ^h, ^^J, Ms, and Aqcd: 

Aqcd 
Qt ' 



1) 


Aqcd 


Ms 


2) 


Ms _^ 

9 ' 






3) 


Aqcd _ 







A, 



QCD 

Q 



(9) 



where L = In(iy), and the series in the exponent makes 



Any Aqcd / A* J power correction can be taken as a cross- 
term between types 1) and 2) for the purpose of enumer- 
ation. The type 1 power corrections are enhanced by the 
presence of the soft scale and are encoded by the moments 
ilk ~ ^QCD tl^*^ function. Type 2 are kinematic 
power corrections that occur because of the expansion 
about small r, and can be computed with perturbation 
theory. The importance of these first two types depends 
on the region considered in Eq. ([5]), with all terms in type 
2 becoming leading order for the far-tail region. Type 3 
are non-enhanced power correction that are of the same 
size in any region. There are also cross-terms between 
the three types. 

In our analysis we keep all power corrections of types 1 
and 2, and the dominant terms of type 3. Our treatment 
of the nonsingular cross section also includes cross-terms 
between 1 and 2 in a manner that is discussed in Sec. HID] 
For the different thrust regions we display the relevant 
terms kept in our analysis in Tab. [lib. The nonsingu- 
lar cross section corrections fully account for the power 
corrections of type 2. The factor [Aqqu / (Qt)]'' in the 
peak region denotes the fact that we sum over all type 1 
power corrections from the leading soft function. In the 
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a) Perturbative Corrections 
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b) Nonperturbative Corrections with Q,k ~ Aqcd 

peak (any k) tail and far-tail (fc = 0, 1, 2) 



das 
dr 

dr 

dr 

p.c. 



i In^r Ofc 
r 

i £ f \ 



Aqcd 



TABLE IL Perturbative and nonperturbative corrections included in our analysis, a) (left panel) Loop orders j for perturbative 
corrections of 0(0?^). Here cusp, non-cusp, and 7^'^ refer to anomalous dimensions, while matching, nonsingular, and the gap 
subtraction S refer to fixed-order series. For convenience in our numerical analysis we use the four-loop beta function for the Oa 
running in all orders displayed, b) (right panel) Nonperturbative corrections included in da/dr with implicit sums over i and 
k. All powers f2fc/(Qr)* can be included in the peak region with the function S*™""*, while only a fixed set of power correction 
parameters are included in the tail and far-tail regions. The row labeled p.c. shows the scaling of the the first power correction 
that is not entirely determined by the earlier rows and hence yield corrections to Eq. (Q. 



tail and multijet regions we only consider the first three 
orders: k=0 (partonic result), k=l (power correction in- 
volving r^i) and , k=2 (power correction involving f22)- 
Here k = 2 terms are used in our error analysis for our 
simultaneous fit to asijnz) and fli. The leading power 
correction that is not fully captured in all regions is of 
type 3, and are of ©(asAqco/Q)- Since our analysis 
is dominated hy Q ~ raz or larger, parametrically this 
gives an uncertainty of 

"^1 ^ ^ ^ 0.3% (10) 
L Us J p.c. Q 

in our final fit (taking Aqcd = 0.3 GeV to obtain the 
number here). This estimate has been validated by run- 
ning our fits in the presence of an additional a.sAQCD/Q 
power correction.^ 

C. Singular Partonic Distribution 

The singular partonic thrust distribution d(Ts/dr con- 
tains the most singular terms oc ai In*^ (r)/r and ai5{T) 
that arise from perturbation theory. Using SCET one 
can derive a factorization theorem for these terms which 
allows for the resummation of the logarithmic terms to 
all orders in perturbation theory. In massless QCD the 
factorization formula for the perturbative corrections in- 
volving Us reads 

— ^ W - Q E "^0 H'q{Q, Uh{Q, m, f^) / ds ds' 



To perform this test we include an as(/ins)Ai/(3 correction in 
the normalized thrust cross section, vary Ai = ±1.0 GeV, and 
perform our default fit to as{mz) and Hi as described in Sec. lVIl 
This variation causes only a ±0.1% change to these fit parame- 
ters, which is smaller than the estimate in Eq. mOI I. 




Here CTq the total partonic e+e~ cross section for 
quark pair production at tree level from a current of type 
/ = {uv, dv, bv, ua, da, ha\ as explained below. Large logs 
are summed by the renormalization group factors Uh be- 
tween the hard scale and /x, Uj between the jet scale and 
/i, and Ug between the soft scale and ji. The choice of /i 
is arbitrary and the dependence on /.i cancels out exactly 
when working at any particular order in the resummed 
expansion. Short distance virtual corrections are con- 
tained in the hard function Hq. The term is the 
thrust jet function. The term 5*^^'' is the partonic soft 
function and the 5(i?, /is)-dependent exponential imple- 
ments the perturbative renormalon subtractions. There 
are four renormalization scales governing the factoriza- 
tion formula, the hard scale ^ Q, the jet scale ^j, 
the soft scale and the renormalon subtraction scale 
R. We have R iis properly sum logarithms related 
to the renormalon subtractions, and there is also a renor- 
malization group evolution in R. The typical values for 
^j, fis, and R depend on r as discussed in Sec. IIIII 

The total tree level partonic e+e" cross section CTq = 
(TQ{Q,mz,Tz) depends on the cm. energy Q, the Z- 
mass, and Z-width, and has six types of components, 
cr^", ctq"", crjj", a^", erg", crg", where the first index de- 
notes flavor, M = up + charm, d = down + strange, 
and b = bottom, and the other index denotes produc- 
tion through the vector (v) and axial-vector (a) cur- 
rents. For QCD corrections we have the hard functions 
= i/^" = H^" = ir§', H^", H^", and i?^", where the 
vector current terms do not depend on the flavor of the 
quark. For massless quark production the axial-vector 
hard functions differ from the vector due to flavor sin- 
glet contributions. All six Cg's and Hq's arc relevant for 
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the implementation of the &-mass and QED corrections. 
Since we use data taken for energies cfose to the Z pole 



we adopt ij [q 



Q Tz/mz) as the Z-boson propa- 



gator which is the form of the width term used for thrust 
data analyses. The modifications of Eq. (|lip required to 
include QED effects are discussed in Sec. Ill HI 

The hard factor Hq contains the hard QCD effects 
that arise from the matching of the two-jet current in 
SCET to fuU QCD. For = Q we have H^{Q,Q) 

1 + hj[as{Q)/4:TrY , and the full hard function with 

\n{^H/Q) dependence is given in Eq. (|A6|) . For the fla- 
vor nonsinglet contributions where the final-state quarks 
are directly produced by the current one can obtain 
the matching coefficient from the on-shell quark vector 
current form factor, which is known to Ofa?) [H, l48l - 
[52}. Converting the bare result in Ref. [53 (see also 
Refs. [11,111]) to the MS scheme and subtracting l/efj^ 
divergences present in SCET graphs, the three-loop non- 
singlet constant, which is one of the new ingredients in 
our analysis, is 
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8998.080, which is the value 



For ny = 5 we have /13 
used for our analysis.'' 

The axial-vector hard functions Hq^ and Hq^ are 
equal to Hq up to additional singlet corrections that 



The analytic expression for /13 in Eq. l|12p is consistent with 
Eq. (7.3) given in Ref. [H. 




wmm 



t,b 



FIG. 2: Two-loop singlet correction to the axial current. Its 
cuts contribute to the hard coefficient and nonsingular terms. 

enter at ©(a^) and 0{al). The fact that the SCET 
hard functions have these singlet corrections was dis- 
cussed in Ref. [55|. At 0(0^) only the axial- vector cur- 
rent gets a singlet correction. It arises from the axial- 
vector anomaly, from suitable cuts of the graph shown in 
Fig. [5] where each axial current is connected to a triangle. 
Summing over the light quarks u,d,s,c gives a vanish- 
ing contribution from this graph, but it does not vanish 
for heavy quarks due to the large bottom-top mass split- 
ting [5^ . Since for the Qs we consider top-pairs are never 
produced, the required terms can be obtained in the limit 
mh/mt — > 0. For the axial current the hard correction 
arises from the bb cut and gives Hq^ = Hq'' = Hq, and 



^singlet ^ where 



hin). (13) 



Here r, 
Ref. [5 



= Q'^/{Am1) and the function l2[rt) from 
is given in Eq. \A7\ . Throughout our analysis 
we use rat = 172 GeV. 77q"^ is a percent level correc- 
tion to the cross section at the Z peak and hence is non- 
negligible at the level of precision of our analysis. (The 
uncertainty in the top mass is numerically irrelevant.) 
At ©(ajj the singlet corrections for vector currents are 
known |42j , but they are numerically tiny. We therefore 
neglect the 0{a'l) vector current singlet corrections to- 
gether with the unknown 0{a'l) singlet corrections for 
the axial- vector current. Likewise we do not account for 
Oia'l) singlet corrections to the nonsingular distributions 
discussed in Sec. Ill El 

The full anomalous dimension of Hq is known at three- 
loops, 0{a^s) [H, Ullll^. It contains the cusp anoma- 
lous dimension, responsible for the resummation of the 
Sudakov double logarithms, and the non-cusp anoma- 
lous dimension. To determine the corresponding hard 
renormalization group factor Uh at the orders N'^LL' 
and N-^LL we need the ©(a^) cusp anomalous dimen- 
sion Tg"*'^ which is still unknown and thus represents a 
source of theory error in our analysis. We estimate the 
size of r™'''' from the order [1/1] Fade approximant in 
Us built from the known lower order coefficients, which is 
within 13% of the two other possible Fade approximants, 
[0/2] and [0/1]. For our theory error analysis we assign 
200% uncertainty to this estimate and hence scan over 
values in the range F™'''' = 1553.06 ± 3016.12. 

The thrust jet function Jt is the convolution of the two 
hemisphere jet functions that describe coUinear radiation 
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in the t and — t directions, 



Jr{s,fJ.)= / ds' J {s' , fi) J {s — s' , fl) 



^ oo 

^— Jn[as{fl)]Cn{s/^l'). (14) 

Here the coefficients J„ are muhiphed by the functions 

■ In" X 



C-i{x) = d{x) 



Cn{x) 



(15) 



where n > 0. Here £„>o(x) are the standard plus- 
functions, see Eq. (|A17|) . At 0{al) only J_i(as) through 
J5(as) are nonzero. The results are summarized in 
Eq. (|A16p . In SCET the inclusive jet function is defined 
as 



J{Qr+,^Ji) 
Im 



AttN^Q 



I d^xe'^-- (0|T{x„(0)^Xn(a^)}|0) 



(16) 



where the x„ are quark fields multiplied by coUinear Wil- 
son lines. The hemisphere jet function has been com- 
puted at 0{as) [H, [531 and O(a^) [HO]. Its anomalous 
dimension is known at three loops, and can be obtained 
from Ref. [6l|. At the order N^LL' we need the 0{al) 
corrections to the jet function. From the anomalous di- 
mension we know the logarithmic terms, Jq to J_5 in 
Eq. at three loops. In the non-logarithmic term J_i 
at 0{al) there is an unknown coefficient (which we de- 
fine as the constant non-logarithmic 3-loop coefficient in 
the position space hemisphere jet function). We estimate 
a range for from the largest value obtained from the 
three Fade approximations for the position space hemi- 
sphere jet function that one can construct from the avail- 
able results. This gives j 3 = ± 3000 for the range of 
variation in our theory error analysis. We note that for 
the 0{a^g) coefficient the corresponding Fade estimate 
/i3 = ± 10000 covers the exact value given in Eq. (fT2)) . 

The renormalization group factors of the thrust jet 
function C/J and thrust soft function Ug sum up large 
logs involving the jet and the soft scales. The required 
cusp and non-cusp anomalous dimensions are fully known 
at three-loops, but again there is dependence on the four- 
loop cusp anomalous dimension r™^P. This dependence 
is included when we scan this parameter as described 
above in our description of the hard evolution. 

The hadronic thrust soft function Sr describes soft ra- 
diation between the two jets. It is defined by 



Sr{k,^i) = -l(o|trrJy„<5(fc. 



id)Y^Yl\0), (17) 



where Yn = Fn(0) and Y^ = Yfl(O) are defined below 
Eq. (|6]). The soft function factorizes into a partonic per- 
turbative part S'p^''* and a nonperturbative part S*™""^, 
Sr = S'P'*'' ® 5'^°'^, as discussed in detail in Sec. HTD] 



This factorization has already been used above in Eqs. ([4]) 
and pT|) . 

At the partonic level the soft function is 

00 

s^\k,^i)^- J2 sMfi)]Cr^{k/fi), (18) 



where S-i to 6*5 arc the only nonzero coefficients at 
0{al), and Cn{x) is defined in Eq. ([T5l) . Results for 
these Skicts) are summarized in Eq. (jA14p . S*?™* was 
calculated at 0{as) in Ref. At 0{a^J the non- 

logarithmic correction was determined in Refs. [lO, 
using numerical output from EVENT2 [H, E^. The 
numerical constant that appears in the non-logarithmic 
0{al) term S-i is referred to as S2 (which is defined as 
the constant 2-loop coefficient in the logarithm of the po- 
sition space soft function). We use S2 = —39.1 ±2.5 [63 |. 
and this uncertainty is taken into account in our theory 
error analysis. * The anomalous dimension of the soft 
function is a linear combination of the anomalous dimen- 
sions of the hard and jet functions which can be obtained 
from the consistency conditions [2^ |4l| . As for the jet 
function we need the 0{al) corrections to S'p^''*. From 
its anomalous dimension we know the logarithmic terms 
at three loops, namely Sq to S5 in Eq. The only 

unknown is the 0{a^) non- logarithmic correction in S-i, 
referred to as S3 (which is defined as the constant non- 
logarithmic term in the logarithm of the position space 
hemisphere soft function). Just like for the constant 
we estimate a value for S3 from the largest value obtained 
from the three possible Fade approximations to the posi- 
tion space soft function that one can construct from the 
available results. This yields the range S3 = ± 500, 
which we scan over in our theory error analysis. 

As already mentioned, in Ref. [l^l an analytic expres- 
sion for the resummed singular thrust distribution was 
presented. Their derivation relies on the Laplace trans- 
form of the jet and soft functions. In our analysis we 
have derived the resummed cross section using two inde- 
pendent procedures, performing all convolutions either 
in momentum space (as presented in App. or in 
Fourier space. These two approaches have been imple- 
mented in two independent codes and we have checked 
that they give exactly the same results. We note that the 
Fourier transform method is equivalent to the Laplace 
procedure used by Bechcr and Schwartz in Ref. (20| 
through a contour deformation, and we find agreement 
with their quoted N^LL formula including matrix ele- 
ments and anomalous dimensions. Furthermore, we also 
agree with their result for the fixed-order singular terms 
up to 0(af). 

In summary, the singular terms in the thrust factor- 
ization theorem are known at N'^LL order, up to the 
unknown constant Fg"^^ . The effect of the cusp anoma- 
lous dimension at 4-loops is much smaller than one might 



Note that in Ref. 62| our S2 was called si. 
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estimate, so for numerical purposes the cross section is 
known at this order. The constants S3 and J3 only en- 
ter for our N'^LL' order. For the singular terms they 
predominantly affect the peak region with spread into 
the tail region only due to RG evolution. Thus in the 
tail region the numerically dominant N'^LL' terms are 
all known. The uncertainties from r™^P. s^, and are 
discussed more explicitly in Sec. IVIII 



D. Oi and Nonperturbative Corrections 

In this section we discuss nonperturbative corrections 
to the thrust distribution included in our analysis, as 
shown in Tab.HIb (right panel). We focus in particular on 
those associated to the first moment parameter Vli. Our 
analysis includes the operator product expansion (OPE) 
for the soft function in the tail region, and combining per- 
turbative and nonperturbative information to smoothly 
connect the peak and tail analyses. We also discuss our 
treatment of nonperturbative corrections in the far-tail 
region, and for the nonsingular terms in the cross sec- 
tion. 

In the tail region where k Qt ^ Aqcd we can per- 
form an operator product expansion of the soft function 
in Eq. At tree level this gives [fillli] 

Sr{k,fi) = S{k) - S'{k) 2ni + ... . (19) 

where the nonperturbative matrix clement fli is defined 
in the MS scheme as 

= F^(o)y„(o)zay„t(o)F;(o)|o) . (20) 

Dimensional analysis indicates that Cti ^ Aqcd- When 
the OPE is performed beyond tree level we must add 
perturbative corrections at a scale ^ ~ A; to Eq. ([TO]). 
The first operator in the OPE is the identity, and its 
Wilson coefficient is the partonic soft function. Thus 
6{k) — > S^'^'-'^{k, fi) when the matching of the leading 
power operator is performed at any fixed order in pertur- 
bation theory. Here we derive the analog for the Wilson 
coefficient of the Cli matrix element and prove that 

Sr{k,pi) = Sr''{k)- ^ ' 2n, + .... (21) 

dk 

This result implies that the leading perturbative correc- 
tions that multiply the power correction are determined 
by the partonic soft function to all orders in perturba- 
tion theory. The proof of Ec^. ((2T|) is given in App. |B] 
The uniqueness of the leading power correction Cli to all 
orders in the perturbative matching can be derived fol- 
lowing Ref. j65| , and we carry out an all orders matching 
computation to demonstrate that the Wilson coefficient 
is dS'P^''(fc)/dfc. At first order in Qi/k «: 1 Eq. (^1]) 
shows that the perturbative corrections in the OPE are 
consistent with a simple shift to Srik — 2fli,fi). This 



type of shift was first observed in the effective coupling 
model (H. 

To smoothly connect the peak and tail regions we use 
a factorized soft function [3, [l^ [s^ 

Sr{k, /i) = Jdk' SP^'\k - k\ y.) 5r'i(fc') , (22) 

where iS'p^'^' is a fixed-order perturbative MS expression 
for the partonic soft function, and 5'™°'^ contains the non- 
perturbative ingredients. In the tail region this expres- 
sion can be expanded for k' <^k and reduces to precisely 
the OPE in Eq. ^ with the identification 

= j dk' k' S'^°'^{k') , (23) 

and normalization condition / dfc' S^'°'^{k') = 1 [ll|. All 
moments of S'^°'^{k') exist so it has an exponential tail, 
whereas the tail for S'P^''*(fc) is a power law. In the 
peak region the full nonperturbative function S'™°'^(fc) 
becomes relevant, and Eq. ([22]) provides a nonperturba- 
tive function whose fi dependence satisfies the MS renor- 
malization group equation for the soft function. In posi- 
tion space the convolution in Eq. (|22[) is a simple product, 
making it obvious that Eq. provides a completely 
general parametrization of the nonperturbative correc- 
tions. The complete basis of functions used to parame- 
terize S'™°'^(fc) in the peak region is discussed in Sec. IIVI 
The expression in Eq. also encodes higher order 
power corrections of type 1 from Eq. ^ through the 
moments 2' f2i = J dk k"^ S^°'^{k), which for tree level 
matching in the OPE can be identified as the matrix 

elements (0|trF^(0)y„(0) (i9)' r„t(o)F!(0)|0)/7Vc. For 
i > 2 perturbative as corrections to the soft function 
OPE would have to be treated in a manner similar to 
App. [B] to determine the proper Wilson coefficients, and 
whether additional operators beyond the powers {iOy 
start contributing. The treatment of perturbative correc- 
tions to these higher order nonperturbative corrections is 
beyond the level required for our analysis. 

Using Eq. the hadronic version of the singular fac- 
torization theorem which involves Sr immediately yields 
Eq. (|lip and the first term in Eq. Q . The conversion of 
^partj-^'j from MS to a renormalon-free scheme is 

discussed in Sec. IIIFI 

Next we turn to the effect of Aqcd power corrections 
on the nonsingular terms in the cross section in Eq. Q. 
The form of these power corrections can be constrained 
by factorization theorems for subleading power correc- 
tions when T ^ 1, and by carrying out an OPE analysis 
for power corrections to the moments of the thrust dis- 
tribution. In the following we consider both of these. 

Based on the similarity of the analysis of power cor- 
rections to thrust with those in i? — )■ Xgj [S^, [68| . 
the factorization theorems for the nonsingular correc- 
tions involves subleading hard functions, jet functions 
and soft functions. They have the generic structure 
i/f '''^ {Q,T,x,)(E) 4"^ {sj ,Xi)(E) si'''' {Qt, Sj/Q) , where the 
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Xi and Sj are various convolution variables. Here si''^ in- 
cludes the leading order soft function in Eq. PT)) as well 
as power suppressed soft functions. Neglecting nonper- 
turbative corrections the nonsingular cross section yields 
terms we refer to as kinematic power corrections of type 
2 in Eq. (|9]). If we do not wish to sum large logs in 
the nonsingular partonic terms, they can be treated in 
fixed-order perturbation theory and determined from the 
full fixed-order computations. In the tail region these r- 
suppressed terms grow and become much more important 
than the Aqcd/Q power corrections of type 3 from sub- 
leading soft functions. In the transition to the far-tail 
region, near r = 1/3, they become just as important as 
the leading perturbative singular terms. In this region 
there are large cancellations between the singular and 
nonsingular terms (shown below in Fig. [7]), and one must 
be careful with the treatment of the nonsingular terms 
not to spoil this. 

We require the nonsingular cross section terms to yield 
perturbative corrections at leading power in Aqcd that 
are consistent with the fixed-order results and with mul- 
tijet thresholds. Our treatment of power corrections in 
the nonsingular terms is done in a manner consistent with 
these goals and with the OPE for the first moment of the 
thrust distribution. To achieve this we use 



moment of the thrust distribution. Eq. (|4]) yields 



dr 



where diTns/dr is the partonic nonsingular cross section 
in fixed-order perturbation theory, whose determination 
we discuss in Sec. Ill El Eq. (p4|) is independent of the 
rcnormalization scale /Xns order by order in its series 
expansion in as(//ns)- The convolution with the same 
S'™°'^(fc') as the singular terms allows the perturbative 
corrections in das/dr + da^s/dr to smoothly recombine 
into the fixed-order result in the far-tail region as re- 
quired by the multijct thresholds. Eq. yields the 
second term in Eq. (U). We will treat the conversion of 
f2i and S"^°'^ to a renormalon-frec scheme in the same 
manner as for the singular cross section, which again for 
consistency requires a perturbative subtraction for the 
partonic da^s/dT that we treat in Sec. IIIFI 

Note that Eq. (|24[) neglects the fact that not all of 
the T dependence in da^s/dT must necessarily be con- 
voluted with 8™°'^. This causes a deviation which is 
^ asTAQCD/{QT) and hence is at the same level as other 
power corrections that we neglect. The largest uncer- 
tainty from our treatment of da^s/dr arises from the 
fact that we do not sum Inr terms, which would require 
anomalous dimensions for the subleading soft and hard 
functions for these nonsingular terms. These logs are 
most important in the peak region, and less relevant in 
the tail region. The size of missing higher order nonsin- 
gular terms such as log enhanced terms will be estimated 
by varying the scale /Xns- 

Our setup is also consistent with the OPE for the first 



dr r- 



dcr 
d7 



dr T 



df 



d(7„ 



df 



2n, 



(25) 



where the ellipses denote C(as Aqcd/Q) and 
0{Aq^Y)/Q^) power corrections. In App. [C] we 
demonstrate that a direct OPE computation for the 
thrust moment also gives the same result, and in 
particular involves precisely the same matrix element 
f2i at this order. The theoretical expression in Eq. 
simultaneously includes the proper matrix elements 
that encode power corrections in the peak region, tail 
region, and for moments of the thrust distribution. 
This implies a similar level of precision for the multijet 
region. Although Eq. (|4]) does not encode all a^AQCD/Q 
corrections, it turns out that the ones it does encode, 
involving numerically give an accurate description of 
the multijet cross section. (This is visible in Fig. [TSIand 
will be discussed further in Sec. IVIII ) This agreement 
provides additional support for our treatment of non- 
perturbative corrections in the nonsingular cross section 
in Eq. 



E. Nonsingular Distribution 

The nonsingular partonic thrust distribution d(Tns/dr 
accounts for contributions in the thrust distribution that 
are kinematically power suppressed. We write 



d(7n 



dr 



f{r/f), (26) 



with the same superscript / notation for different cur- 
rents as in Eq. pT|) . The presence of the d{R,jis)- 
dependent exponent arises because 5™°'^ depends on fii 
and we use the same renormalon-free definition for fli as 
for the singular terms. In our numerical evaluation we 
integrate by parts so that the d/dr derivative acts on 
^mod 1^ This exponent is discussed in detail in 

SecEFl 

In this section we discuss our determination of the 
functions in pure QCD with massless quarks, while 
the generalization to include mi, effects is discussed in 
Sec. Ill Gl and to include QED effects in Sec. Ill HI For pure 
QCD there is one function f^^^^ = /^^ = ^ = /qcd ^r 
the vector current, and functions f^^^ = /q°j, and /^^j 
for the axial-vector currents. In general is the par- 
tonic fixed-order distribution where the singular terms 
which are already contained in diTs/dr are subtracted to 
avoid double counting. Setting the renormalization scale 
fJ-ns = Q they have the form 



/qcd(T, 1) 



■./2(r) 



0(^,l)=/qcd(^,l): 



(2^) = 
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FIG. 3: 0{as) nonsingular thrust distribution. 
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FIG. 4: ©(a?) 

nonsineular thrust distribution. 
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FIG. 5: 0{ai) nonsingular thrust distribution. For simplicity 
we only show the data binned with 0.01 bin size. 



/nc"d(^,l)=./ncd(r,l) + 



: /singlet (t ,rt) : 



(27) 



where here as = as{Q) and rj = Q'^/{Am^). The re- 
quired results for /^(r, Uns/Q) can be obtained by shift- 
ing UsiQ) to as(/ins) using the fixed-order relation be- 
tween these couplings at 0{al). 



The full 0{as) partonic thrust distribution has been 
known analytically for a long time [6^ . For the one- loop 
nonsingular distribution it gives 



3t(t- 1) 



6 + 6 r - 4) log ( i - 2 



3t 



[3 + 41og(T)]. 

(28) 



This result is plotted in Fig. |3l The kink at r — 1/3 
appears because the full one-loop distribution vanishes 
at this value with a nonzero slope, and there is an exact 
cancellation between the fixed-order singular and non- 
singular one-loop expressions. For r > 1/3 the one-loop 
nonsingular distribution is precisely the negative of the 
one-loop fixed-order singular distribution. 

The 0{a'^) and 0{al) QCD distributions are available 
in numeric form from the Fortran programs EVENT2 [gI, 
m and EERAD3 Q (see also Ref. |H^), respectively. 
These programs are used to derive results for our /2(t) 
and fsir) nonsingular distributions in a manner dis- 
cussed below. At ©(ofg) there is also the singlet correc- 
tion fsingietiT,r) for the axial- vector contribution aris- 
ing from the large bottom-top mass splitting. The 
three-parton quark-antiquark-gluon cut from Fig. [5] con- 
tributes to the nonsingular distribution, and we have in- 
cluded this contribution analytically [70|. The formula 
for /singlot (r, r) is given in Eq. ()A30|) . There is also a 
contribution from the four-parton cut. Its contribution 
to .fsinsict(T,r) is unknown, but it is tiny for the total 
cross section [5g and can therefore be safely neglected. 

At 0{al) we use linear binned EVENT2 results for 
T > 0.095 and log-binning results for t < 0.095 each 
obtained from runs with 10"'^° events and infrared cut- 
off yo = 10"*. For r > 0.095 (using a 0.005 bin size) 
the resulting statistical uncertainties in the nonsingular 
distribution are always below the percent level and negli- 
gible and we can use an interpolation of numerical tables 
for I2{t). For r < 0.095 the singular terms dominate the 
distribution which leads to large cancellations and an en- 
hancement of the statistical uncertainties. Here we use 
the ansatz /2(t) = Y^i=o ^^i' t + t Yl^i=2 ^^"^ 
the coefficients and hi to the EVENT2 output, includ- 
ing the constraint that the integral over the full distribu- 
tion reproduces the known 0{a1) coefficient for the total 
cross section. The result has the form /2(t) -I- £2^2 ('''), 
where /2 represents the best fit and 82 is the 1-sigma error 
function with all correlations included. The term e2 is a 
parameter which we vary during our as-^i fit procedure 
to account for the error. Here /2 and 82 also depend on 
the coefficient S2 in the partonic soft function St which is 
known only numerically. In Fig. |4]we plot the EVENT2 
data we used, along with our /2(t) with S2 = —39.1. The 
dashed curves show the result for £2 = ±1, with the re- 
gion inbetween corresponding to the 1-sigma error band. 

For the determination of /s at 0{a^g) we implement a 
similar approach as for /2 , using results from EERAD3 Q 
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computed with 6 x 10^ events for the three leading color 
structures and 10'' events for the three subleading ones, 
using an infrared cutoff yo = 10^'"'. We employ linearly 
binned results with 0.01 bin size for r > 0.315 (keep- 
ing the statistical error below the percent level) and with 
0.005 bin size for r < 0.315. For the fit for r < 0.315 our 
ansatz function has the form /3(r) = '^^^iCihi^'T and 
the result has the form /s (r) + 63 ^3 (r) , with being the 
best fit and ^3 the 1-sigma error function. The constant 
£3 is the analog of £2 and is varied in the error analysis. 
We note that and S3 depend on the constant S2 and 
on the constants S3 and J3 that account for the unknown 
non-logarithmic terms in the 0{al) soft and jet func- 
tions. This dependence is included in our error analysis. 
In Fig. M we plot the EERAD3 data with bin size 0.01, 
along with our hir) with S2 = -39.1, /13 = 8998.08, 
j3 = S3 = 0. The dashed curves show the result for 
£3 = ±1, with the region inbetween corresponding to the 
1-sigma error band. 

In our analysis we use the values —1,0, 1 for £2 and £3 
to account for the numerical uncertainties of our fit func- 
tions in the small r region. The nonsingular partonic dis- 
tribution depends on one common renormalization scale 
/ins which is varied in our theory error analysis as given 
in Sec. ITTIl 



F. Gap Formalism 

The partonic soft function S^'^'^^{k) computed pertur- 
batively in MS has an O(Aqcd) renormalon ambigu- 
ity. The same renormalon is present in the partonic MS 
thrust cross section with or without resummation. This 
is associated with the fact that the partonic threshold at 
fc = in S'P'*'^'(fc) is not the same as the physical hadronic 
threshold for the distribution of soft radiation that occurs 
in S'r(fc). One can see this explicitly in the large- /3o ap- 
proximation, where it is associated to a pole at u = 1/2 
in the Borel transform [T6| 



B 



16C 



-5/6 Q 

(29) 



This result shows that S'P'^'''(fc) in the MS scheme suf- 
fers from the renormalon ambiguity for all fc > 0. The 
MS matrix element Jli defined in Eq. ((20|) also has an 
©(Aqcd) renormalon ambiguity. Together, the renor- 
malon in this power correction and in the perturbative 
series for S^^'^*'{k) combine to give a soft function Sr{k) 
that is free from this O(Aqcd) renormalon. If left un- 
subtracted this renormalon ambiguity leads to numerical 
instabilities in perturbative results for the thrust distri- 
bution and in the large order dependence for the deter- 
mination of the soft nonperturbative function 5'™°'^. In 
this section we resolve this problem by switching to a new 
scheme for Qi. This scheme change induces subtractions 



on dcrf"^''*/dr that render it free of this renormalon. We 
start by reviewing results from Ref. [l6j . 

Consider a class of soft nonperturbative functions with 
a gap parameter A, which only have support for k > 
A, so S'^°'i(fc) ^ S^°'^{k - 2A). Here the MS moment 
relation in Eq. becomes 



2A+ /dfcfcS'™°'^(fc) = 2^1 



(30) 



where A accounts for the complete renormalon ambigu- 
ity contained in f2i. We can now obtain a renormalon- 
free definition for i^i by splitting A into a nonpertur- 
bative component A{R,fis) that is free of the ©(Aqcd) 
renormalon, and a suitably defined perturbative series 
6(R, fis) that has the same renormalon ambiguity as fli. 
The parameter A is scheme and renormalization group 
invariant, while A and 6 individually depend on the sub- 
traction scale R and in general also on the soft scale us- 
Writing 



A ^ AiR, fis) + S{R, ^ls) , 



(31) 



the factorization of perturbative and nonperturbative 
components in Eq. ((22|) becomes 



Sr{k,fis) =jdk' S^'''\k-k'-2S,ns) S^°'^{k' -2A) 
' dk' [e-2'5*S'P^^t(fc-fc',/is)j5;"°'^(fc'-2A). (32) 



Here the exponential operator induces perturbative sub- 
tractions (in powers of as{ns)) on the MS series in 
^partj-^-j ^Y^g^^ render it free of the renormalon. This expo- 
nential modifies perturbative results for the cross section 
in the manner we have shown earlier in Eqs. (|11|) and 
The convolution of the nonsingular cross-section 
with 8^°'^ in Eq. ([24]) now becomes 



dk' 



da„ 



(33) 



Furthermore, with Eq. ([5^ the result in Eq. ([50)) becomes 



2AiR,^ls)+ JdkkS^'°''ik)^2ni{R,ns), 



(34) 



where here fli{R,fis) is renormalon-free. Combining 
Eqs. (|34p and (pO)) we see that the scheme conversion 
formula from MS to the new scheme is 



ni{R, ^is) = ^1 - S{R, ^ls) ■ 



(35) 



Thus, the precise scheme for fli{R, ns) is specified by the 
choice of the subtraction series J(i?,/is)- Note that in 
general the gap parameter A is an additional nonpertur- 
bative parameter that can be determined together with 
other parameters in the function 5'™°'* from fits to exper- 
imental data. However, in the tail region the power cor- 
rections are dominated by a single parameter, fli{R, fJ-s), 
which encodes the dependence on A. 
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In Ref. [G^l a convenient scheme for 5{Rjj^) was de- 
rived (based on a scheme proposed in Ref. [7l|) where 



d \n{ix) 



In Sr(x, ^) 



x = (iB£'<E )- 



(36) 



Here Sr{x, fi) is the position space partonic soft function, 
and the fact that we write this result for S^- rather than 
for the hemisphere soft function explains the extra factor 
of 1/2 relative to the formula in Ref. [62|. The cutoff pa- 
rameter R, having mass dimension 1, is a scale associated 
with the removal of the infrared renormalon. To achieve 
the proper cancellation of the renormalon in Eq. p2p one 
has to expand 6{R,iJ,s) together with S^^'^*'{k, ^s) order 
by order in a^ips)- The perturbative series for the sub- 
traction is 



SiR,fis) = e'^ 



(37) 
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FIG. 6: The running of Q,i{R,R) with R = i?(r), plotted as 
a function of r for Q = 35, 91.2, 207 GeV. 



where the 5i>2 depend on both the adjoint Casmir Ca 
3 and the number of light flavors in combinations that 
are unrelated to the QCD beta function. For five light 
flavors the one, two, and three-loop coefhcients are 



7^ = 1615.42228 + 54.6195541 sa , 



(40) 



5i{R,iis) = -0.848826Lii, 

hiR, IJ-s) = -0.156279 - 0.46663Li? - 

S3{R,fis) = 0.0756831 + 0.01545386 S2 



0.517864L|j, 



0.777219L1 



0.421261L|, 



0.622467Lij 

(38) 



with Lj^ = \n{fis/R)- We will refer to the scheme defined 
by Eq. ((36|) as the R-gap scheme for fli. 

From the power counting Qi ^ Aqcd one expects that 
a cutoff R^ 1 GeV should be used, such that fii ~ Aqcd 
and perturbation theory in as{R) remains applicable. 
We refer to this as the power counting criterion for R. 
Since in the tail region ^5 ~ Qt ^ 1 GcV the factors 
of L]^ in Eq. (f551) arc then large logs. To avoid large 
logarithms in the subtractions Si{R,^s) it is essential to 
choose R ^ fis, so that the subtraction scale R is de- 
pendent on T much like the soft scale fis- We refer to 
this as the large-log criterion for R. To resolve the con- 
flict between these two criteria, and sum the large logs 
while keeping A{R,jiq ~ R) renormalon- free, we make 
use of R-cvolution [391 . l40j . Formulas for the gap case 
were given in Ref. [62| and are reviewed here. In this 
scheme A(i?, /i) satisfies an R-RGE and /i-RGE 



dR 
d 



n+l 



dfj, 



(39) 



so that 7^ = — 2e'''^r'^"''P[Q;s]. For five fiavors the anoma- 
lous dimension coefficients up to three loops are 



70^=0, 7f 



while the coefficients F^'^p are given in Eq. (|A26|) . The 
solution of Eq. ^ at N'^LL is 

A{R, ^) = A(i?A, A^a) + i?eT^c^[F™^P, fi, R] 
+ i?Ae'^^w[r™^P,i?A,MA] 
+ A[^l^D^''^[a,{R),as{RA)], (41) 



where the resummed ^[F'^^^p, /i, ^q] is given in Eq. (|A23[) 
and the resummed D^'"'^ [as(i?), as(i?A)] is given in 
Eq. (jA3ip . Both the gap subtraction and R-evolution 
equations at C(a^) depend on the constant S2 which 
we vary within its errors in our theory error scan. In 
our analysis, when quoting numerical results, we always 
use the parameter A(i?A,A'A) at the reference scales 
^A ~ ^J'A ~ 2 GeV to satisfy the power counting cri- 
terion for R. We then use Eq. (|4ip to run up to the 
scale i? /i5 in order to satisfy the large-log criterion. 
The precise R value is a function of t, i? = R{t), and 
given in Sec. IIIII with our discussion of the profile func- 
tions. The RGE solution for A(i?, ^5) in Eq. (gT]) yields a 
similar solution for a running r2i(i?, /ig) using Eq. (|34p . 
In Fig. ini we show the result for the running ili{R,R) 
with the boundary value ^i{Ra, (J'a) = 0.323 GeV. The 
anomalous dimension and R(t) profile function cause an 
increase in the size of the power correction for increasing 
r and for increasing Q. 

Note that our R-gap subtraction scheme differs from 
the subtractions in the low-scale effective coupling model 
of Ref. [35|, which is not based on the factorization of 



-43.954260 , 



the soft large angle radiation but on the assumption that 
the ©(Aqcd) renormalon ambiguity is related entirely 
to the low-energy behavior of the strong coupling a^. In 
the effective coupling model the subtractions involve log- 
arithms, hi(fi/fii), where fi is the usual renormalization 
scale of perturbation theory and fij is the low- momentum 
subtraction scale, which is set to /i/ = 2 GeV. The scale 
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FIG. 7: Components of the pure QCD cross section. Here 
fii = 0.35 GeV and as(mz) = 0.114. 



fij plays a role very similar to the scale R in the R- 
gap scheme. These logarithms are the analogs of Ln in 
Eq. ([55)1 and, since n (x Q these logarithms also become 
large. In the effective coupling model an appropriate re- 
summation formalism for large logs in the subtractions 
remains an open question. 

In Fig. [7] we plot the absolute value of four components 
of our cross section for our complete QCD result at N'^LL' 
order in the R-gap scheme at Q = mz- The cross sec- 
tion components include the singular terms (solid blue), 
nonsingular terms (dashed blue), and separately the con- 
tributions from terms that involve the subtraction coef- 
ficients Si, for both singular subtractions (solid red) and 
nonsingular subtractions (dashed red) . The sum of these 
four components gives the total cross section (solid black 
line). The subtraction components are a small part of 
the cross section in the tail region, but have an impact 
at the level of precision obtained in our computation. In 
the peak region at very small r the solid red singular 
subtraction grows to be the same size as the solid blue 
singular term, and is responsible for yielding a smooth 
positive definite total cross section. In both the peak 
and tail regions the singular cross section dominates over 
the nonsingular cross section. But as we approach the 
threshold r ~ 1/3 for the far-tail region they appear with 
opposite signs and largely cancel. This is clear from the 
figure where individually the singular and nonsingular 
lines are larger than the total cross section in this region. 
The same cancellation occurs for the singular subtraction 
and nonsingular subtraction terms. 



G. Bottom Mass Effects 

In this work we implement bottom mass effects 
using the SCET factorization framework for massive 
quarks [13, EH • We include mh-depcndence in the kine- 



matics, which starts at tree level, and in the 0{as) cor- 
rections in the partonic singular and nonsingular distri- 
butions. We also account for the resummation of large 
logs and for hadronization effects in the mf,-dependcnt 
terms. The mass dependent factorization theorem im- 
plies that the rcnormalization group summation of loga- 
rithms is identical to the one for massless quarks, and 
that all power corrections of type 1 from Eq. ^ are 
described by the nonperturbative soft function S"^°^ al- 
ready defined for the massless case [l3, EH. We have 
already indicated this with the convolution Adat/dr (8) 
^mod gi^own in Eq. Since for the numerical analysis 
in this work we fit to data in the tail region, where Qt > 
6 GcV; and since the massive quark thrust factorization 
theorem implies for the soft scale fis ^ Qt > 6 GeV, 
we do not have to account for any flavor threshold in 
the rcnormalization group evolution and can always use 
71/ = 5. The mass dependent factorization theorem fur- 
ther implies that the only nontrivial TOf,-dependence in 
the singular distribution arises in the thrust jet function. 
Thus the jet scale ^ Qy/r ^ uib for the region of our 
fit and we use the MS bottom mass rni,{fj:j) to parame- 
terize the mi, corrections with rhi,{mb) = 4.2 GeV as our 
input value. Using the MS mass rather than the pole 
mass avoids the appearance of large higher order effects 
related to the ©(Aqcd) polc mass renormalon. 

We implement the partonic bottom mass corrections 
as an additive term to the massless partonic N'^LL' cross 
section. These corrections come from the production of 
bottom quarks by the virtual 7 or Z, 

Ad<Tb ^dab_ da^''=° ^^^^ 
dr dr dr ' 

where both d(Tf,/dr and d(7™''~"/dr are computed at 
NNLL. Because the effect of rhb 7^ in Ad<T;,/dr is ex- 
pected to be a percent level correction to the tail cross 
section, we anticipate that the NNLL level of precision 
suffices. (This is also justified a posteriori by the rela- 
tively small effect of the nib corrections on our fit results.) 

An important aspect in the discussion of the finite 
quark mass effects is in which way hadron and heavy 
quark masses need to be accounted for in the definition 
of thrust in Eq. ([T]). In the experimental analyses Monte 
Carlo generators are used to convert the actual measure- 
ments to the momentum variables needed to compute r, 
and this conversion depends on hadron masses. Since the 
final state stable hadrons are light, these effects are re- 
lated to nonperturbative physics. Theoretically they are 
therefore implicitly encoded within our fit of the nonper- 
turbative corrections. In the partonic theoretical com- 
putation light hadron masses are neglected in the com- 
putation of the r distribution, and it is consistent to set 
J2i \Pi\ = Q i'^ the denominator of Eq. ([1]). 

To understand how the heavy quark masses affect the 
definition of thrust in Eq. ([1]) we recall that the partonic 
computation relies on the inclusive nature of the mea- 
surements and that, experimentally, only light and long- 
lived hadrons reach the detectors and are accounted for 
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in the Pi momenta that enter in computing r. Thus for 
heavy hadrons containing bottom (or charm) quarks, it 
is their hght and long-hved hadronic decay products that 
enter the particle sum Due to energy conservation 
it is therefore necessary to set \pi\ = Q in the denom- 
inator of the thrust definition of Eq. ([T]) for the leading 
power partonic computations involving heavy quarks. On 
the other hand, due to three-momentum conservation, it 
is consistent to use the heavy quark three-momentum in 
the numerator of Eq. ([T]) for the partonic computations. 
This makes the partonic thrust computations involving 
heavy quarks simple because we do not need to explic- 
itly account for the heavy quark decay in the calculations. 
Together with the relation \pi\ = Q in the denomina- 
tor of Eq. (ITJ this induces a shift of the observed thrust 
distribution for b quarks to larger r values. Comparing 
to the massless quark situation, the small-r cndpoint is 
moved from to 



rr" = l-\/l-4m2/g 



(43) 



where here rhb = fh{,{lJ-j)- At tree level this shifts 
6{t) 6{t — t"""). For the fixed-order result at 
0{as) the three-jet cndpoint is moved from 1/3 to 
r^3jet = 5/3 _ 4/3v'l -3?Ti^/Q2. At leading order in 
fhl/Q^ < 1 we have r'™" = 2fhl/Q^ + 0{fn^/Q^) 
and r^^j''* = 1/3 + 2ml/Q^ + 0{rhl/Q^), so the shift 
is the same for both endpoints. Numerically, for mf, = 
4.2 GeV and Q = (35, 91.2, 207) GeV, r is shifted by 
(0.029,0.004,0.0008). This shift is also observed experi- 
mentally in fiavor tagged thrust analyses (t^ - It^ ]. 

In the following we outline the method used to com- 
pute the partonic diTh/dr. Like for the massless case 
the distribution is divided into singular and nonsingular 
parts 



dob ^ d^l^do^ 
dr dr dr 



(44) 



The implementation of the bottom mass effects into the 
singular distribution do-^/dr follows the NLL' analysis in 
Ref. [Ill, except that the evolution in the present work 
is incorporated fully at NNLL order and that the exact 
partonic threshold at r = r^"'" is accounted for, 

X Jrb{s', rht, fij) Uj{s - s', ^, /ij) j dk Ug{k, /i, fis) 



X e ^ « a 



(MS-polc mass scheme change terms) , 



(45) 



where crg(a;) = crg^Vl- 4^2(1 + 2x^) + cr[;"(l - Ax^)^/"^. 
Perturbative bottom mass effects in the soft function 
start at two loops, so at 0{as) S'P'^''* remains unchanged. 
Since we have fhf,/Q 1, only the thrust jet function 
for bottom quark production, JTb{s,fhi,, fi) [75|, receives 



modifications from the finite rrif,. These modifications 
lead to a shift of the partonic threshold of the thrust jet 



function from invariant mass p — to p 



In 



Jrb{s,fhb, ^jl) the variable s = p^ — m^, and the presence 
of the mass leads to r™'" in Eq. (|45| . It also gives a 
more complicated form for 0{as) corrections in Jrb in- 
volving regular functions of m^/ s in addition to singular 
terms oc 5{s) and ^ti' [s / fi^) / [s / familiar from the 
massless quark jet function. More details and explicit 
formulae can be found in Refs. [T^. l4l|. 

The bottom quark mass effects in the nonsingular par- 
tonic distribution da^^/dr are more complicated since 
finite mass effects at 0{as) differ for vector and axial- 
vector current induced jet production, 



da," 



rUb Mns 



) 



Tn-b ^J.ns 

Q ' Q 



Q ' Q 

+ (MS-pole mass scheme change terms) 



(46) 



In our analysis we implement analytic expressions for the 
nonsingular functions and The full 0{as) dis- 
tributions for T > can be obtained from integrating 
the known double differential bb energy distribution for 
vector-induced and axial- vector-induced production, re- 
spectively, see e.g. Refs. [zl, [t^].^ The corresponding 
0{as) coefficient of the S{t — r™'") term is obtained us- 
ing the one-loop correction to the total bb cross section as 
a constraint. To determine the nonsingular distributions 
/^'" we proceed much like for the massless case and sub- 
tract the singular contributions expanded to C(as) from 
the full 0{as) distribution. Further details and explicit 
formulas for /^"'" will be given in a future publication. 



H. QED Corrections 

For the electroweak corrections to the thrust distri- 
bution we can distinguish purely weak contributions and 
QED effects. The dominant effects to jet production from 
the purely weak interactions are given by virtual one-loop 
corrections to the hard Wilson coefficient Hq. Since the 
contribution of the singular thrust distribution dag/dr 
dominates in the t ranges we use for our fits as well as in 
the total cross section crtot = / dr dcr/dr (see Fig. [7|) , the 
purely weak corrections largely drop out when the distri- 
bution is normalized to the total cross section. This is 
consistent with the explicit computations carried out in 
Refs. [isl . |43 | , where purely weak corrections were found 
to be tiny. In our analysis we therefore neglect purely 
weak effects. 



® Results for bottom mass corrections at O(a^) were determined 
in Refs. [78l480l |. but are not used in our analysis due to the small 
effect the bottom mass corrections have in our fits. 
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For QED corrections the situation is more comphcated 
because, apart from virtual effects which again largely 
cancel in the normalized distribution, one also has cor- 
rections due to initial state and final state radiation. In 
addition, one has to account for the fact that the treat- 
ment of QED effects in the thrust measurements depends 
on the experiment. In general, using Monte Carlo sim- 
ulations, all experimental data were corrected to elim- 
inate the effects from initial state radiation. However, 
they differ concerning the treatment of final state pho- 
ton corrections, which were either eliminated or included 
in the corrected data sets. In Sec. IVII we review infor- 
mation on the approach followed by the various experi- 
mental collaborations. Since many experiments did not 
remove final state radiation, we have configured a ver- 
sion of our code that adds final state photons and QED 
Sudakov effects, and does so on an experiment by exper- 
iment basis. A parametric estimate of the potential im- 
pact of these QED effects on the measurement of as{mz) 
is ^ -0.2AAaem/iCFas) ^ -1%, where 0.244 is the av- 
erage of the square of the electromagnetic charges for the 
five lightest fiavors. 

We implement the leading set of QED corrections to all 
components that go into the main factorization formula 
of Eq. Q in the massless quark limit counting ^ ctg 
to make a correspondence with Tab.HIland remembering 
to include cross terms such as terms of 0{acmCts)- Ex- 
ceptions where QED corrections are not included are the 
gap subtraction S{R,fis) and the R-evolution equation 
for the gap parameter A. This is because QED effects do 
not lead to ©(Aqcd) infrared renormalon ambiguities. 
Most of the required QED results can be obtained in a 
straightforward manner from modifications of the known 
QCD corrections. 

Our implementation of QED effects is briefly described 
as follows: For the evolution of the strong coupling we in- 
cluded the ©(a^Qiow) corrections to the QCD beta func- 
tion. There are also effects from the evolution of the 
QED coupling aami^J■) which we define in the MS scheme. 
In the beta function for the QED coupling acm we ac- 
count for the dominant 0{a'^^^) and the next-to- leading 
0{a^j^as) contributions. For the full singular partonic 
distribution which includes both QCD and QED effects 
we have 

d(T f 
' = Q^a'„ HQ{Q,HH)UH{Q,fiH,fJ-) dsds' 

I •' 



dr 



(47) 



s 

Q 



where all factors now depend on the index / due to their 
dependence on the electromagnetic charges q^="^."'' = 
-f2/3 and qi=dv,da,hvM ^ implement one- 

loop QED corrections in the hard factor _ffg, the jet 
function and the soft functions S^^^^^ . In the renor- 
malization group evolution factors t/^, ?7j^, C/^^ we ac- 
count for the one-loop QED corrections to the cusp and 



the non-cusp anomalous dimensions. In the nonsingu- 
lar partonic distribution dtT,is/dr the same approach is 
employed. Here the C'(aom) contributions that are anal- 
ogous to the 0[pLs) terms are included by writing the full 
functions to be used in Eq. ([26]) as 



I\2 



8tt 



-Mr). (48) 



The 1% parametric estimate and the moderate size of 
the QED effects we observe from the results of our fits 
justifies the neglect of higher order QED effects. A more 
precise treatment of QED effects is also not warranted 
given the level of accuracy of the Monte Carlo generators 
used to correct the experimental data. More details and 
explicit formulae for the QED corrections discussed here 
will be given in a future publication. 



III. PROFILE FUNCTIONS 

The factorization formula for the singular partonic dis- 
tribution d(Ts/dr in Eq. pT|) is governed by three renor- 
malization scales, the hard scale ^h, the jet scale /ij, 
and the soft scale fis- To avoid large logarithms appear- 
ing in the corrections to the hard coefficient Hq , the jet 
function Jr and the soft function Sr, the corresponding 
scales must satisfy the following theoretical constraints 
in the three r regions: 



1) peak: IJ-h ^ Q , ^ v^AqcdQ , 

2) tail: ^^h Q , fJ-j QVt' , 

3) far-tail: fin = f-tj = l^i-s ^ Q ■ 



Ms ^ AqcD , 
(49) 



In the peak region, where the full nonperturbative func- 
tion 5'™°'^ is relevant we have fin ^ fij ^ IJ'S Aqcd- 
In the tail region, where the nonperturbative effects are 
described by a series of moments of the soft function we 
have fiH fJ-j ^ fJ-s ^ Aqcd- To achieve an accurate 
theoretical description, we resum logarithms of r in the 
peak and tail region where fin, fJ-j^ and fis are separated. 
Finally, in the far-tail region the partonic contributions 
are described by usual fixed-order perturbation theory, 
and a proper treatment of fixed order multijet thresholds 
requires that the three n parameters merge close together 
in the far-tail region and become equal at r = 0.5, with 
^J'H = = Ms ^ <3 ^ Aqcd- Thus in the far-tail region 
logarithms of r are not summed. The merging of hh, 
fij^ and fis in the far-tail region is of key importance for 
the cancellations between singular and nonsingular cross 
sections shown in Fig. [T] To obtain a continuous fac- 
torization formula that is applicable in all three regions 
we use T-dependent renormalization scales, which we call 
profile functions. These are smooth functions of r which 
satisfy the theoretical constraints listed in Eq. (j^Hl)- 

In addition to the three renormalization scales of the 
singular partonic distribution there are two more scales, 
/ins and R. The renormalization scale /ins governs the 
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pcrturbative series for the function contained in the 
nonsingular partonic distribution d(T„s/dr. The subtrac- 
tion scale R arises when we implement the gap sub- 
tractions in the R-gap scheme for Qi that remove the 
C(Aqcd) renormalon contained in the MS soft func- 
tion. This R also corresponds to the endpoint of the 
R-cvolution for A{R, fis) given in Eq. (pij) . To avoid 
large logarithms in the subtraction S{R, fJ.s)y the value of 
R needs to be chosen of order fis and is therefore also a 
function of r. 

The factorization formula Q is formally invariant un- 
der 0(1) changes of the profile function scales, that is, 
changes that do not modify the hierarchies. The residual 
dependence on the choice of profile functions constitutes 
one part of the theoretical uncertainties and provides 
a method to estimate higher order perturbative correc- 
tions. Wc adopt a set of six parameters that can be var- 
ied in our theory error analysis which encode this residual 
freedom while still satisfying the constraints in Eq. ([1^)) . 

For the profile function at the hard scale, we adopt 



100 r 




0.4 ^ 0.5 



FIG. 8: Profile functions for the renormalization scales ^j{t), 
/is(T), and subtraction scale R{t) that appear in the factor- 
ization theorem. Shown are results for the central parameter 
values aX Q = mz- 



(50) 



where en is a free parameter which we vary from 1/2 to 
2 in our theory error analysis. 

For the soft profile function we use the form 



-0.848826 ln(/is/i?) is nonzero (see Eq. 
fore use the form 



). We there- 



i?(r) 



^0 + l-llT + I^L2t'^, 

A*s(t), 



< r < ii, 
ti < r < 0.5. 



(53) 



Ms(t) = ^ br + d, 



< r < ii, 

tl<T< t2, (51) 
i2 < T < i. 



Here, ti and t2 represent the borders between the peak, 
tail and far-tail regions. ^0 is the value of at t = 0. 
Since the thrust value where the peak region ends and 
the tail region begins is Q dependent, we define the Q- 
independent parameter ni = ti{Q/lGeW). To ensure 
that (r) is a smooth function, the quadratic and linear 
forms are joined by demanding continuity of the function 
and its first derivative at r = <i and t = t2, which fixes 
& = 2 {fiH - m)/{t2 - ti + ^) and d = [^0(^2 + 5) - 
I^H ti] / (^2 ~ ^1 + ^) • In our theory error analysis wc vary 
the free parameters ni, t2 and /xq. 

The profile function for the jet scale is determined by 
the natural relation between the hard, jet, and soft scales 



t^j{r) = 1 



(52) 



The term involving the free C'(l)-parameter ej imple- 
ments a modification to this relation and vanishes in the 
multijet region where r ~ 1/2. We use a variation of ej 
to include the effect of such modifications in our estima- 
tion of the theoretical uncertainties. 

For the subtraction scale R the choice R = /is(T) en- 
sures that we avoid large logarithms in the Si{R,iis) 
subtractions for the soft function. In the peak re- 
gion, however, it is convenient to deviate from this 
choice so that the 0{as) subtraction term 5i{R, ^s) ~ 



Imposing continuity of R{t) and its first derivative at 
T ^ ti requires fii = {2d — 2i?o + bti)/ti and ^2 = 
{~d + i?o)/tf . The only free parameter is Rq which sets 
the value of i? at t = 0. We take Rq = 0.85/io to give the 
one loop subtraction Si{R, fis) the appropriate sign to 
cancel the renormalon in the peak region. Since our focus 
here is not the peak region, we leave further discussion 
of the appropriate choice of Rq to a future publication. 

In our theory error analysis we vary to account 
for our ignorance on the resummation of logarithms of 
T in the nonsingular corrections. We account for the 
possibilities 



Mns(T) 




Ms(t")), 



ns = 1, 
ns = 0, 
= — 1. 



(54) 



We do not include the choice fins = fJ-s since we find that 
the choice of this small scale enhances the nonsingular 
contributions in an unnatural way. 

In total, we have introduced six free parameters which 
we vary to account for renormalization scale uncertain- 
ties. In our analysis we use the following central val- 
ues and variations: /_to = 2 



0.25 



-1-0.05 
^0.05i 



0^ 



en 



:o.5 GeV, m : 
= 2'' with h 



'+3 



t2 = 



Ot\ and 



Us = (—1, 0, 1). In Fig. [5] we show the form of the profile 
functions for Q = mz = 91.2 GeV and all profile param- 
eters at their central values. The dashed lines represent 
the functions Q\/t and Qt which were the central choices 
for /ij(T) and ijls{t) used in Ref. [1^, but which do not 
meet in the multijet region. In order for our profile for 
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/is(T) to join smoothly onto fin and /^./(t) it is necessary 
for /Z5(t) to have a slope ~ 2Qt in the tail region. Since 
In 2 is not large our profiles sum the same Inr's as with 
the choice in Ref. [201, but satisfy the criteria necessary 
to treat the multijet thresholds. ^° 



IV. NONPERTURBATIVE MODEL FUNCTION 

The soft nonperturbative function S™°'^{k) parameter- 
izes the dominant nonperturbative hadronic effects in the 
thrust distribution. It describes the hadronization con- 
tributions that arise from how soft hadrons that are ra- 
diated in between the jets enter the thrust variable in 
Eq. ([1]). It is normalized, has the property Sr{0) = 0, 
is positive definite and has support for fc > 0. To keep 
the representation of 5™°'* as much as possible indepen- 
dent of a particular analytic parametrization we adopt 
the approach of Ref. [s^ and write the soft nonperturba- 
tive function as a linear combination of an infinite set of 
basis functions which can in principle describe any func- 
tion with the properties mentioned above. The model 
function we use has the form 



for our purposes, so we set ci = 0. For this case 



5r'(fc,A,{cO) = ^ 



TV 



(55) 



where the basis functions are [38 



fniz) 



2z3(2n + 1) 



g[z) = ^ (3 - e-'^^ (3 + 12z + + 322'"^)) - 1, (56) 

and Pn are Legendre polynomials. For c| = 1 the 
norm of 5'™°'^(/c) is unity, VIq — 1. The choice of basis in 
Eqs. (j55|) and (j56p depends on specifying one dimension- 
ful parameter A which is characteristic of the width of 
the soft function. With N = 00 the parameter A would 
be redundant, but in practice we truncate the sum in 
Eq. ([55]) at a finite N, and then A is effectively an addi- 
tional parameter of the model function. 

In this work we fit to experimental thrust data in 
the tail region where the predominant effects of the 
soft model function are described by its first moment 
ni(A, A, {ci}). As explained below, we use the second 
moment n2{X,S,{ci}) to validate our error analysis and 
confirm the validity of neglecting this parameter in the 
fit. Since in the tail region the exact form of the soft 
model function is not relevant, we take N ^ 2 setting 
c„>2 = 0. Variations of the parameter ci are highly cor- 
related with variations of A and are hence not necessary 



In Ref. [23! where NLL resummation is achieved by exponentia- 
tion, the log resummation is turned off at a predefined threshold 
Tmax with the log-R method [2ll | . In this approach the transition 
to fixed order results in the multijet region differs from ours. 



A 



0.201354coC2 + 1.10031c^] 



n2=A^ + AX[cl + 0.201354coC2 -I- 1.10031c^] 
A2 



■ [l.25co + 1.03621coC2 + 1.78859( 



2J 7 



(57) 



and the normalization condition + = 1 can be 
used to eliminate cq > 0. Recall that in the soft 
model function in the factorization theorem we must use 
Sf°'^{k-2A[R,^ls),\{c^}) where R = R{t) and fis = 
fJ-siT) are determined by the profile functions. When we 
quote numbers for parameters we use A ~ A(i?A,MA) 
and hence f2i.2 = ^^i,2(^A, /^a) with reference scales 
/iA = ^A = 2 GeV. The running between the scales 
{R,fis) and (i?A,MA) is determined by Eq. (|4T|) . 

For our default fit in the tail region only the parameter 
ill is numerically relevant so without loss of generality we 
can take cq = 1, C2 = 0, and set A(i?A,MA) = 0.05 GeV. 
In this case all higher moments fin>i are determined as 
a function of J7i and A. For example we have f22 ~ 
(A^ — 2Ar2i -t- for the second moment. 

In Sec. IVIII we analyze the dependence of our fit re- 
sults on changes of il2- Because C2 has a rather strong 
correlation to il2, we implement these variations by 
using Eq. (|57p and setting C2 to nonzero values. In this 
case we can hold fli fixed by a suitable choice of A for a 
given C2- 

To obtain results from our code that do not include 
nonperturbative corrections we can simply turn them off 
by setting S'^'°'^{k) = S{k) and A = (5 = 0. 



V. NORMALIZATION AND CONVERGENCE 

The experimental data is normalized to the total num- 
ber of events. In our prediction we therefore need to 
normalize the distribution to the total cross section, i.e. 
we have to calculate {1 / a)da / dr . Since the factorization 
formula in Eq. Q is valid for all thrust values we have 
the option to use either the integral of our da/dr distri- 
bution for the norm, or the available fixed-order result 
for the total hadronic cross section. 

The fixed-order total cross section is 



Had 



i?"" = i?Had + Ra+it^ I{rt) , i?"" = Rmd+Rv 
37r^ 



(58) 



Here i?Had is the pure QCD cross section for mass- 
less quarks, Ra.v are mass corrections depending on 
rrih/Q, and I{rt) is the isosinglet correction from the axial 
anomaly and large top-bottom mass splitting [HBl- Set- 
ting ji = Q the QCD cross sections for massless quarks 
at three loops is 



^H, 



ad 



1 -f- 0.3183099 as{Q) + 0.1427849 al{Q) 
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- 0.411757a^(Q) . (59) 

We refer to the review in Ref. [8l| for a discussion of the 
fixed-order hadronic cross section. We note that the 
series for the fixed-order hadronic cross section exhibits 
an excehent and fast convergence. At 0{al) the pertur- 
bative uncertainty is much below the permille level and 
hence entirely negligible for the purpose of our analysis. 

In the R-gap scheme in pure QCD, from a numerical 
analysis at Q — niz, we find at N^LL' order that the 
integrated norm of the thrust distribution for the default 
setting of all theory parameters (see Tab. lIIip gives about 
0.99 

""tot ^i'^s)- However we also find that the pertur- 
bative uncertainty of the integrated norm (determined by 
the theory scan as described in Sec. IVip is about ±2.5%, 
which is substantially larger than for the fixed-order cross 
section. This larger uncertainty is due to the perturba- 
tive errors of the thrust distribution in the peak region. 
At N^LL' order we therefore employ the fixed order cross 
section to normalize the thrust distribution we use for the 
fits. 

At the lower orders in the R-gap scheme (N'^LL, 
NNLL', NNLL, NLL') we find that the integrated norm 
for central theory parameters is more appropriate since 
the order-by-order convergence to crf^ is substantially 
slower than that of the rapid converging fixed-order QCD 
result in Eq. ([59]) . Again we find that the large perturba- 
tive uncertainties in the peak region render the perturba- 
tive errors of the integrated norm larger than those of the 
fixed-order norm. We therefore evaluate the integrated 
norms at the lower orders with the theory parameters 
fixed at their default values (see Tab. This means 

that to estimate the theoretical errors in our fits to exper- 
imental data at orders below N^LL' in the R-gap scheme, 
we vary the theory parameters only for the distribution 
and not for the norm computation. In the MS scheme for 
fii we also adopt the integrated norm at all orders. When 
we evaluate the thrust distribution with log-resummation 
but without nonperturbative effects we use the same nor- 
malization choices as for the R-gap scheme, which makes 
comparison to earlier work in Sec. lIXI easier. For the situ- 
ation where the cross-section is evaluated at fixed-order, 
without resummation or nonperturbative effects, we use 
the appropriate fixed order normalization at each order. 

As discussed in Sec. |Vll to compare with the binned 
experimental data we integrate our theoretical expression 
for the distribution (l/(T)(d(T/dr) over each bin [Ti,r2]. 
A potential alternative is to use theoretical results for the 
cumulant 

I](r)= rdr'-^(r'). (60) 

Here one sums large logs of r rather than r', and 
the SCET based cumulant has r dependent profiles, 
S(t, /^^(t)). The presence of fJ-i{T) implies that the 
derivative of the cumulant is not precisely equal to the 
distribution, 

±^r,^,,{r))^i^{r,t,,{r)) (61) 



The difference coming from the second term in Eq. (|6ip 
can be numerically important for certain observables. To 
test this we consider using for the cross-section integrated 
over the bin [Ti,r2] the theoretical expression 

S(r2,^,(f2)) - i;(ri,^,(fi)) , (62) 

and will examine several choices for fi,2. 

One simple possibility is to use fi — ti and f2 = T2, 
so that S(t2, //i(T2)) — E(ri,^i(ri)) is used. In this case 
there is a spurious contribution from outside the [ti,T2] 
bin associated to the second term in Eq. (pT|) . 

I](ri,^j(T2)) - S(ri,^i(Ti)) (63) 

where the ~ holds under the approximation that the 
derivative do not change very much across the bin. With 
our default setup the deviation of this simple choice 
for the cumulants from our integrated result for the 
distribution is 2% to 8% for r G [0.1,0.3], bin-size 
T2 - Ti = 0.01, and Q = 91.2 GeV." In the far-tail 
region ti S [0.3,0.45], where the cross-section becomes 
small, the deviation grows from 8% to 1000%. These 
deviations are dominated by the spurious contribution. 
The size of the spurious contribution is not reduced by 
increasing the bin-size to T2 — ri = 0.05, and is only 
mildly dependent on Q. Any choice in Eq. (j62p where 
fi ^ T2 leads to a spurious contribution from r € [0,ri]. 

If we instead use ti = f2 = {ti+T2)/2 then the spurious 
contribution is identically zero. In this case the difference 
between Eq. ((62)) and our integrated thrust distribution is 
reduced to 0.5% for n e [0.1,0.3] and for n G [0.3,0.45] 
grows from 0.5% to only 20%. Although dramatically re- 
duced, the difference to the integrated distribution in the 
far-tail region is still quite sizeable. This discrepancy oc- 
curs because only for the distribution (1 / a){(la / dr) can 
the (r) profile functions be constructed such that they 
satisfy exactly the criteria discussed in Sec. IIIII Due to 
the above issues, and since the binned datasets are in- 
tended as representations of the thrust distribution, we 
have determined that our approach of integrating the 
thrust distribution is conceptually the best. 

In the rest of this section we discuss the perturbative 
behavior of the thrust distribution in the tail region. The 
values of the physical parameters used in our numerical 
analysis are collected in Eq. (|A4[) . For our lower order 
fits we always use the four-loop beta function in the run- 
ning of the strong coupling constant, as mentioned in the 



For the profile functions used by Becher and Schwartz [20II . dis- 
cussed in section [iXl this deviation has similar size but opposite 
sign. 
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FIG. 9: Theory scan for errors in pure QCD with massless quarks. The panels are a) fixed-order, b) resummation with no 
nonperturbative function, c) resummation with a nonperturbative function using the MS scheme for Cli without renormalon 
subtraction, d) resummation with a nonperturbative function using the R-gap scheme for Qi with renormalon subtraction. 



caption of Tab. HT] Furthermore, we always consider five 
active flavors in the running and do not implement bot- 
tom threshold corrections, since our lowest scale in the 
profile functions (the soft scale /ig) is never smaller than 
6 GcV in the tail where we perform our fit. 

In Fig. [9] we display the normalized thrust distribution 
in the tail thrust range 0.15 < r < 0.30 at the differ- 
ent orders taking as{mz) = 0.114 and ^^i{Ra, (J-a) ~ 
0.35 GeV as reference values, and neglecting rub and QED 
corrections. We display the case Q = mz where the 
experimental measurements from LEP-I have the small- 
est statistical uncertainties. The qualitative behavior of 
the results agrees with other cm. energies. The colored 
bands represent the theoretical errors of the predictions 
at the respective orders, which have been determined by 
the scan method described in Sec. I VII 

In Fig. Oi we show the 0{as) (hght /yellow), O(a^) 



(medium/purple) and 0{al) (dark/red) fixed-order 
thrust distributions without summation of large loga- 
rithms. The common renormalization scale is chosen 
to be the hard scale fin- In the fixed-order results the 
higher order corrections are quite large and our error es- 
timation obviously underestimates the theoretical uncer- 
tainty of the fixed-order predictions. This panel including 
the error bands is very similar to the analogous figures 
in Refs. Q and @. This emphasizes the importance of 
summing large logarithms. 

In Fig. IHb the fully resummed thrust distributions at 
NLL' (yellow), NNLL (green), NNLL' (purple), N^LL 
(blue) and N'^LL' (red) order are shown, but without 
implementing the soft nonperturbative function 5'™°'^ or 
the renormalon subtractions related to the R-gap scheme. 
The yellow NLL' error band is mostly covered by the 
green NNLL order band, and similarly the purple NNLL' 
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band is covered by the blue N^LL one. Moreover the blue 
N'^LL band is within the purple NNLL band. Compared 
to the fixed-order results, the improvement coming from 
the systematic summation of large logarithms is obvious. 
In particular we see that our way of estimating theoret- 
ical uncertainties is appropriate once the logarithms are 
properly summed. At N-^LL and at N-^LL' order the rela- 
tive uncertainties of these resummed thrust distributions 
in the tail region r e [0.1,0.3] are about ±7.8% and 
±4.6%, respectively. 

The results shown in Fig. |9l: are very similar to panel b 
but now include also the soft nonperturbative function 
^mod -^yitiiout renormalon subtractions, where f2i is de- 
fined in the MS scheme. In the tail region the soft non- 
perturbative function leads to a horizontal shift of the 
distribution towards larger thrust values by an amount 
St (X 2fli/Q. This is clearly visible by comparing the 
values at r = 0.15 where the curves intersect the j/-axis. 
Concerning the uncertainty bands and the behavior of 
predictions at the difi^erent orders the results are very 
similar to those in panel b. 

Finally, in Fig. I^li we show the results with summa- 
tion of large logarithms including the soft model function 
with renormalon subtractions, where fii is defined in the 
R-gap scheme. In the R-gap scheme the convergence of 
perturbation theory is improved, and correspondingly the 
size of the uncertainties from the same variation of the 
theory parameters is decreased. The decrease of the un- 
certainties is clearly visible comparing the blue N^LL and 
the red N'^LL' uncertainty bands with panel c. The rela- 
tive uncertainties of the thrust distribution at N'^LL and 
at N'^LL' order in the tail region r € [0.1,0.3] are now 
about ±3.4% and ±1.7%, respectively. This improve- 
ment illustrates the numerical impact of the ©(Aqcd) 
renormalon contained in the partonic soft function and 
shows the importance of eliminating the ©(Aqcd) renor- 
malon. 



VI. EXPERIMENTAL DATA AND FIT 
PROCEDURE 

Experimental data for thrust are available for various 
cm. energies Q between 14 and 207 GeV. In our analy- 
sis we fit the factorization formula (|H) in the tail region 
to extract as and fti. As our default data set we use 
the thrust range 6/Q < t < 0.33, and we only employ 
data from Q > 35 GeV. The lower boundary 6/Q removes 
data in the peak where higher order moments become im- 
portant, while the upper boundary of 0.33 removes data 
in the far-tail region where the OsAqcd/Q power cor- 
rections become more important. We take Q > 35 GeV 
since a more sophisticated treatment of b quark effects 
is required at lower energies. The data we use are 
from TASSO with Q = {35,44} GeV [13, AMY with 
Q = 55.2 GeV [H, JADE with Q = {35,44} GeV H, 
SLC with Q = 91.2 GeV [H, L3 with Q = {41.4, 55.3, 
65.4, 75.7, 82.3, 85.1, 91.2, 130.1, 136.1, 161.3, 172.3, 



182.8, 188.6, 194.4, 200.0, 206.2} GeV [73, IH], DEL- 
PHI with Q = {45, 66, 76, 89.5, 91.2, 93, 133, 161, 
172, 183, 189, 192, 196, 200, 202, 205, 207} GeV [13- 
lai, OPAL with Q = {91, 133, 161, 172, 177, 183, 189, 
197} GeV and ALEPH with Q = {91.2, 133, 161, 

172, 183, 189, 200, 206} GeV [9|. (For TASSO and 
AMY we have separated statistical and systematic errors 
using information from the experimental papers.) All 
data is given in binned form, and we therefore integrate 
Eq. (HI over the same set of bins to obtain appropriate 
theory results for the fit to the experimental numbers. 
For the case that either r = 6/Q or r = 0.33 are lo- 
cated within an experimental bin, that bin is excluded 
from the data set if more than half of it lies outside the 
chosen interval. For the Q > mz data we removed five 
bins with downward fluctuations that were incompatible 
at the > 10-sigma level with the cross section implied by 
neighboring data points and other experimental data in 
the same region. The list of these bins is: L3 (136.1 GeV): 
[0.25,0.275], DELPHI (161 GeV): [0.32,0.40], DELPHI 
(183 GeV): [0.08,0.09], DELPHI (196 GeV): [0.16,0.18], 
ALEPH (200 GeV): [0.16, 0.20]. Our defauh global 
data set contains a total of 487 bins. In the numerical 
analysis performed in Sec. IVIII we also examine alterna- 
tive global data sets with different r-ranges. 

The data sets were corrected by the experiments to 
eliminate the QED effects from initial state radiation us- 
ing bin-by-bin correction factors determined from Monte 
Carlo simulations. The primary aim of these corrections 
was to eliminate the effective reduction of the cm. energy 
available for the production of the hadronic final state. 
In addition, in the data sets from the TASSO, L3 and 
ALEPH collaborations the effects from final state radi- 
ation of photons were eliminated, while they have been 
fully included in the data sets from the AMY, JADE, 
SLC, DELPHI and OPAL collaborations. It should also 
be noted that the approaches used by the experiments 
to treat photon radiation were dependent on the cm. 
energy Q. For the Q = mz data any radiation of ini- 
tial state photons is naturally suppressed as the effective 
cm. energy for the hadronic final state gets shifted away 
from the Z pole. Therefore no specific photon cuts were 
applied for the Q ~ mz data prior to the application of 
the bin-by-bin correction factors. For the data taken off 
the Z pole for either Q < mz or Q > mz the effects of 
initial state radiation are substantial and explicit hard 
photon cuts were applied in the data taking prior to the 
application of the bin-by-bin correction procedure. We 
therefore consider the Q = mz data sets as more reliable 
concerning the treatment of QED effects. 

Since the size of the QED effects we find in the mea- 



Four out of these bins lie in our r G [6/Q, 0.33] default fit range. 
If they are included in the default dataset then for our final fit in 
Eq. (|68]l the = 439 increases by +81 and the central fit values 
show a slight decrease to as{mz) = 0.1132 and a slight increase 
to Qi = 0.336 GeV. 
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parameter 


default value 


range of values 


Mo 


2GeV 


1.5 to 2.5 GeV 


ni 


5 


2 to 8 


t2 


0.25 


0.20 to 0.30 


ej 





-1,0,1 


en 


1 


0.5 to 2.0 


Us 





-1,0,1 


S2 


-39.1 


-36.6 to -41.6 


-pcusp 
^ 3 


1553.06 


-1553.06 to H-4569.18 


33 





-3000 to +3000 


S3 





-500 to +500 


£2 





-1,0,1 


£3 





-1,0,1 



35 GeV 



TABLE III: Theory parameters relevant for estimating the 
theory uncertainty, their default values and range of values 
used for the theory scan during the fit procedure. 



suremcnts of as and the soft function moment fii is com- 
parable to the experimental uncertainties (see the results 
and discussions in Sec. IVII[) . a less Monte Carlo depen- 
dent treatment of QED radiation would be certainly war- 
ranted. (See Ref. Q for a recent discussion of QED 
radiation based on full one-loop matrix elements.) How- 
ever, given that the impact of QED corrections we find 
for as and fli is still smaller than the current theoretical 
uncertainties from QCD, we use for our default numeri- 
cal analysis the theory code with QED effects switched 
on, as described in Sec. Ill Hi In Sec. IVHI we also present 
results when QED corrections are neglected for all data 
sets, and for the case when they are neglected only for 
the TASSO, L3 and ALEPH data sets. 

For the fitting procedure we use a x^-analysis, where 
we combine the statistical and the systematic experimen- 
tal errors into the correlation matrix. We treat the sta- 
tistical errors of all bins as independent. The systematic 
errors of the bins are correlated, but - unfortunately - 
practically no information on the correlation is given in 
the experimental publications. Wc therefore have to rely 
on a correlation model. For our analysis we assume as the 
default that within one thrust data set, i.e. for the set of 
thrust bins obtained by one experiment at one Q value, 
the systematic experimental errors are correlated in the 
minimal overlap model used by the LEP QCD working 
group [91I, |93|. In the minimal overlap model the off- 
diagonal entries of the experimental covariance matrix 
for the bins i and j within one data set are equal to 
[min(A^'"', A^-*"')]^, where A^-*^'' are the systematic errors 
of the bins i and j. This model implies a positive correla- 
tion of systematic uncertainties within each thrust data 
set. As a cross check that our default correlation model 
does not introduce a strong bias we also carry out fits 
were the experimental systematic errors are assumed to 
be uncorrelated. Details are given in Sec. I VIII 

To estimate the theoretical errors in the a^-fii plane at 
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FIG. 10: Difference between default cross section and the 
cross section varying only one parameter as a function of 
r. We vary as{mz) by ±0.001 (solid red curves), 2^1 by 
+ 0.1 (dashed blue curves) and C2 by ±0.5 (dash dotted green 
curves). The plot is shown for three different values of the cen- 
ter of mass energy: (a) Q = 35 GeV, (b) Q = 91.2 GeV, (c) 
Q = 206 GeV. 



any order and for any approximation used for the factor- 
ization formula (|4]), we carry out independent fits for 500 
different sets of theory parameters which are randomly 
chosen in the ranges discussed in the previous sections 
and summarized in Tab. IIIIl We take the area covered 
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by the points of the best fits in the Ug-^i plane as the 
theory uncertainty treated hke 1-sigma.^'^ We emphasize 
that this method to estimate theoretical errors is more 
conservative than the error band method [26j employed 
for example in Refs. [lO, [l^]- However, our method re- 
quired considerably more computer power and it was nec- 
essary to use the Tier- 2 centers at Garching and MIT, as 
well as clusters at the MPI and the University of Arizona. 
In Sec. IVIII we also present the outcome of other ways to 
estimate the theoretical error. 

It is an important element of our analysis that we carry 
out global fits to the data from all values of Q > 35 (and 
all experiments) . This is motivated by the strong degen- 
eracy between and fli in the tail region which can only 
be lifted when data from different Q values are simulta- 
neously included in the fits.^"* In Fig. [TU] the difference 
dcr/dr — (dcr/dr)dofauit is displayed for 0.08 < r < 0.30 
and Q = 35, 91.2 and 206 GcV. Here (dcr/dr)dcfauit is 
the cross section for the default setting of the theory pa- 
rameters with as{mz) = 0.114 and fli = 0.35 GeV and 
for da/dr we vary either as{mz) by ±0.001 (sohd red 
curves) or 2f2i by ±0.1 GeV (dashed blue curves) from 
their default values. The figures show that in the tail 
region changes in as can be compensated by changes in 
r^i. This degeneracy makes it impossible to determine 
as and fli simultaneously with small uncertainties from 
tail fits that use data from one Q value (or from a nar- 
row range of Q values). On the other hand, we see that 
the correlation is Q dependent when considering a large 
enough range of Q values. In our fits it is particularly 
important to include, apart from the data from Q = mz, 
the low-energy data from JADE, TASSO, and AMY, and 
the high energy data from the LEP-II experiments. Al- 
though the errors in these analyses are larger than from 
the high-statistics Q = mz run at LEP-I these data sets 
are essential for breaking the degeneracy and simultane- 
ously extracting as and 17 1. 



VII. NUMERICAL ANALYSIS 

Having explained all ingredients of the factorization 
formula ([4]) and the fit procedure we are now in the posi- 
tion to discuss the numerical results of our analysis based 
on a global fit to the experimental data for Q > 35 GeV in 
the tail region. In the tail region the dominant power cor- 
rections are encoded in the first moment fii, see Eq. ([6]), 
so we can determine as{mz) and Qi from a simultane- 
ous fit. In this section we examine in detail the numerical 
results of our fits concerning the treatment of the pertur- 



This corresponds to a 1-sigma error (68% CL) in as as well as 
in Qi. 

The presence of this degeneracy is presumably also related to 
why Monte Carlos that are tuned to LEP data tend to have 
smaller hadronization corrections at Q = mz than at larger Q 
values. See Sec. IIXI 



order a,(mz) (with fif^) a^(mz) (with fif'^^'P) 



NLL' 


0.1203 ±0.0079 


0.1191 ±0.0089 


NNLL 


0.1222 ±0.0097 


0.1192 ±0.0060 


NNLL' 


0.1161 ±0.0038 


0.1143 ±0.0022 


N^LL 


0.1165 ±0.0046 


0.1143 ±0.0022 


N^LL' (full) 


0.1146 ±0.0021 


0.1135 ±0.0009 


N'^LL'(QCD+,Ti6) 


0.1153 ±0.0022 


0.1141 ±0.0009 


N'^LL'(pupo QCD) 


0.1152 ±0.0021 


0.1140 ±0.0008 



TABLE IV: Theory errors from the parameter scan and cen- 
tral values for Qa(mz) at various orders. The N^LL' value 
above the horizontal line is our final scan result, while the 
N'^LL' values below the horizontal line show the effect of leav- 
ing out the QED corrections, and leaving out both the 6-mass 
and QED respectively. The central values are the average of 
the maximal and minimal values reached from the scan. 



order 


Qi (MS) 


f2i (R-gap) 


NLL' 


0.264 ±0.213 


0.293 ± 0.203 


NNLL 


0.256 ±0.197 


0.276 ±0.155 


NNLL' 


0.283 ± 0.097 


0.316 ±0.072 


N^LL 


0.274 ± 0.098 


0.313 ±0.071 


N^'LL' (full) 


0.252 ± 0.069 


0.323 ±0.045 


N^LL'cQCD + nif,) 


0.238 ± 0.070 


0.310 ±0.049 


N^LL'(puro QCD) 


0.254 ± 0.070 


0.332 ± 0.045 



TABLE V: Theory errors from the parameter scan and cen- 
tral values for f^i defined at the reference scales Ra = fJ,A = 
2 GeV in units of GeV at various orders. The N^LL' value 
above the horizontal line is our final scan result, while the 
N'^LL' values below the horizontal line show the effect of leav- 
ing out the QED corrections, and leaving out both the 6-mass 
and QED respectively. The central values are the average of 
the maximal and minimal values reached from the scan. 

bative, hadronization and experimental errors, QED and 
bottom mass corrections and their dependence on the 
choice of the data set. We note that the values quoted 
for r^i in the R-gap scheme are given for reference scales 
i?A = MA = 2 GeV, see Sec. IHFI 

Theory Scan 

In Fig. [TT] the best fit points of the theory parame- 
ters scan in the as-2^li plane are displayed at NLL' 
(brown), NNLL (magenta), NNLL' (green), N^LL (blue) 
and N-^LL' (red) order. The fit results at N'^LL' order 
include bottom mass and QED corrections. In Fig. [TTk 
the results in the R-gap scheme with renormalon sub- 
tractions are shown, and in Fig. Illb the results in the 
MS scheme without gap subtractions are given. 

At each order 500 fits were carried out with the the- 
ory parameters randomly chosen in the ranges given in 
Tab. mil As described in Sec. IVI[ we take the size of the 
area in the as-2D,i plane covered by the best fit points 
as a measure for the theoretical uncertainties. To vi- 
sualize the theoretical uncertainties we have colored the 
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FIG. 11: Distribution of best fit points in the Qs(mz)-2ni and as{mz)-2rii planes. Panel (a) shows results including pertur- 
bation theory, resummation of the logs, the soft nonperturbative function and Qi defined in the R-gap scheme with renormalon 
subtractions. Panel (b) shows the results as in panel a, but with Qi defined in the MS scheme without renormalon subtractions. 
In both panels the respective total (experimental+theoretical) 39% CL standard error ellipses are displayed (thick dark red 
lines), which correspond to 1-sigma (68% CL) for either one-dimensional projection. 
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FIG. 12: Distribution of best fit points in the as(mz)-x^/dof plane. Panel (a) shows the x^/dof values of the points given in 
Fig. Illb . Panel (b) shows the x^/dof values of the points given in Fig. Illb . 



respective areas according to the orders. The fit results 
clearly show a substantial reduction of the theoretical 
uncertainties with increasing orders. Explicit numerical 
results for the respective central values (determined by 
the mean of the respective maximal and minimal values) 



and the theory errors (determined by half of the differ- 
ence between maximal and minimal values) for as and 
r^i are given in Tabs. IIVI and |Vl respectively. We will 
consider these theory errors as 1-sigma. At N'^LL' or- 
der with rii in the R-gap scheme the theory error for 
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FIG. 13: Thrust distribution at N^LL' order and Q = mz 
including QED and mt corrections using the best fit values 
for as(mz) and Q,i in the R-gap scheme given in Eq. (|68[) . The 
pink band represents the perturbative error determined from 
the scan method described in Sec. IVII Data from DELPHI, 
ALEPH, OPAL, L3, and SLD are also shown. 



asjmz) is ±0.0009 compared to ±0.0021 with f^i in the 
MS scheme. Also at NNLL' and N^LL we see that the 
removal of the ©(Aqcd) rcnormalon leads to a reduction 
of the theoretical uncertainties by about a factor of two 
in comparison to the results with Cli in the MS scheme 
without rcnormalon subtraction. The proper treatment 
of the rcnormalon subtraction is thus a substantial part 
of a high-precision analysis for 17i as well as for a^. 

It is instructive to analyze the minimal values for 
the best fit points shown in Fig. [TTJ In Fig. [T^] the dis- 
tributions of the best fits in the ag-xLin/'^'^^ plane are 
shown using the color scheme of Fig. [11] Figure [T2k dis- 
plays the results in R-gap scheme, and Fig. [T2b the ones 
in the MS scheme. For both schemes we find that the 



values and the size of the covered area in the as 



x4in/dof plane systematically decrease with increasing 
order. While the analysis in the MS scheme for fii leads 
to x^jjjj/dof values around unity and thus an adequate 
description of the entire global data set at N'^LL' order, 
we see that accounting for the rcnormalon subtraction in 
the R-gap scheme leads to a substantially improved the- 
oretical description having Xmin/dof values below unity 
aheady at NNLL' and N^LL orders, with the N^LL' or- 
der result slightly lower at x^^j^^/dof ~ 0.91. This demon- 
strates the excellent description of the experimental data 
contained in our global data set. It also validates the 
smaller theoretical uncertainties we obtain for as and fii 
at N^LL' order in the R-gap scheme. 

As an illustration of the accuracy of the fit, in Fig. [T3] 
we show the theory thrust distributions at Q = mz for 
the full N^LL' order with the R-gap scheme for f^i, for 
the default theory parameters and the corresponding best 
fit values shown in bold in Tabs. IIVI and [Vj The pink 



Band 


Band 


Our scan 


method 1 


method 2 


method 


U.UUU4 


O.UOUO 


o.oooy 


0.0016 


0.0019 


0.0021 


0.0018 


0.0021 


0.0034 


0.0018 


0.0026 


0.0046 



N^LL' with nfs'P 
N^LL' with nf^ 
N^LL' without 3^°"^ 
0{al) fixed-order 



TABLE VI: Theoretical uncertainties for as{mz) obtained at 
N^LL' order from two versions of the error band method, and 
from our theory scan method. The uncertainties in the R-gap 
scheme (first line) include renormalon subtractions, while the 
ones in the MS scheme (second line) do not and are therefore 
larger. The same uncertainties are obtained in the analysis 
without nonperturbative function (third line). Larger uncer- 
tainties are obtained from a pure 0(al) fixed-order analysis 
(lowest line). Our theory scan method is more conservative 
than the error band method. 



band displays the theoretical uncertainty from the scan 
method. The fit result is shown in comparison with data 
from DELPHI, ALEPH, OPAL, L3, and SLD, and agrees 
very well. (Note that the theory values displayed are 
actually binned according to the ALEPH data set and 
then joined by a smooth interpolation.) 

Band Method 

It is useful to compare our scan method to determine the 
perturbative errors with the error band metho d 12611 that 
was employed in the analyses of Refs. [20|.[2^|25|. In the 
error band method first each theory parameter is varied 
separately in the respective ranges specified in Tab. Hill 
while the rest are kept fixed at their default values. The 
resulting envelope of all these separate variations with 
the fit parameters as{mz) and Cli held at their best fit 
values determines the error bands for the thrust distri- 
bution at the different Q values. Then, the perturbative 
error is determined by varying as{mz) keeping all the- 
ory parameters to their default values and the value of 
the moment Qi to its best fit value. The resulting per- 
turbative errors of as{mz) for our full N^LL' analysis in 
the R-gap scheme are given in the first line of Tab. IVII 
In the second line the corresponding errors for as{mz) 
in the MS scheme for fii are displayed. The left column 
gives the error when the band method is applied such 
that the as(mz) variation leads to curves strictly inside 
the error bands for all Q values. For this method it turns 
out that the band for the highest Q value is the most 
restrictive and sets the size of the error. The resulting 
error for the N'^LL' analysis in the R-gap scheme is more 
than a factor of two smaller than the error obtained from 
our theory scan method, which is shown in the right col- 
umn. Since the high Q data has a much lower statistical 
weight than the data from Q = mz, we do not consider 
this method to be sufficiently conservative and conclude 
that it should not be used. The middle column gives the 
perturbative error when the band method is applied such 
that the as{mz) variation minimizes a function which 
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FIG. 14: Distribution of best fit points at N^LL' order witii 
Sli in tfie R-gap scheme in pure QCD (light green), including 
7Jit effects (medium blue) and including mi, effects and QED 
corrections (dark red) . Solid circles indicate the central points 
for these three cases. The hollow circle represents the central 
point from the global fit with QED corrections neglected for 
the data from TASSO, L3 and ALEPH, but included for all 
other data sets. 



puts equal weight to all Q and thrust values. This sec- 
ond band method is more conservative, and for the N'^LL' 
analyses in the R-gap and the MS schemes the resulting 
errors are only 10% smaller than in the scan method that 
we have adopted. The advantage of the scan method we 
use is that the fit takes into account theory uncertainties 
including correlations. 

Effects of QED and the bottom mass 

Given the high-precision we can achieve at N^^LL' or- 
der in the R-gap scheme for fii, it is a useful exercise 
to examine also the numerical impact of the corrections 
arising from the nonzero bottom quark mass and the 
QED corrections. In Fig.[T3]the distributions of the best 
fit points in the as-2ili plane at N'^LL' in the R-gap 
scheme is displayed for pure massless QCD (light green 
points), including the bottom mass corrections (medium 
blue points) and the bottom mass as well as the QED 
corrections (dark red points). The distribution of the 
best fit points with bottom mass and QED corrections 
(dark red points) was already shown in Fig. [TTk . The 
large black dots represent the corresponding central val- 
ues. The corresponding numerical results are shown at 
the bottom of Tabs. IV] and El 

We see that the QED and bottom quark mass effects 
are somewhat smaller than the theoretical errors of the 
N'^LL' analysis but not negligible. Moreover wc find that 



the qualitative impact of the QED and the bottom quark 
mass effects is quite intuitive: The nonzero bottom quark 
mass primarily causes a horizontal shift of the thrust 
distribution towards larger r values, since the small-r 
threshold for massive quark production is moved to a 
finite T value. Here this is compensated primarily by a 
reduced value of fii. Concerning QED effects, they cause 
an effective increase of the coupling strength in the final 
state interactions leading primarily to a decrease of as in 
the fit. 

As explained in Sec. I VII the experimental correction 
procedures applied to the AMY, JADE, SLC, DELPHI 
and OPAL data sets were designed to eliminate initial 
state photon radiation, while those of the TASSO, L3 
and ALEPH collaborations eliminated initial and final 
state radiation. It is straightforward to test for the effect 
of these differences in the fits by using our theory code 
with QED effects turned on or off depending on the data 
set. Since our procedure treats data from different 
experiments as uncorrelated it is also easy to implement 
this technically. Using our N^LL' order code in the R-gap 
scheme wc obtain the central values as{mz) — 0.1136 
and fii = 0.318 GcV, indicated by the hollow circle 
in Fig. 1141 Comparing to our default results given in 
Tabs. |IV] and El which are based on the theory code 
were QED effects are included for all data sets, we see 
that the central value for as is larger by 0.0001 and the 
one for Oi is smaller by 0.006 GeV. This shift is substan- 
tially smaller than our perturbative error, and justifies 
our choice to use the theory code with QED effects in- 
cluded as the default code for our analysis. 

Hadronization and Experimental Error 

An important element in the construction of the func- 
tion used for our fit procedure is the correlation model 
for the systematic uncertainties given for the experimen- 
tal thrust bins. The results discussed above rely on the 
minimal overlap model for the systematic experimental 
errors explained in Sec. IVIl The 1-sigma ellipse based 
on the central values of Eq. (|M1) and centered around 
(a^,217i) = (0.1135,0.647 GeV) is shown in Fig. HSlby 
the red solid ellipse. This ellipse yields the experimen- 
tal errors and hadronization uncertainty related to fii 
in our analysis. We find that the size and correlation 
coefficients of the 1-sigma error ellipses at N'^LL' order 
of all fits made in our theory scan are very similar, and 
hence we can treat the theory error and these hadroniza- 
tion/experimental errors as independent. 

The correlation matrix of the red solid error ellipses is 

(i, j = as,2fJi) 



(64) 



3.29(16) • 10-^ -2.30(12) • lO"'"^ GeV ^ 
-2.30(12) • 10-5 GcV 1.90(18) • lO^^ GcV^ ^ 
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FIG. 15: Experimental 1-sigma standard error ellipse (red 
solid) in the as-2r2i plane. The larger ellipse shows the total 
uncertainty including theory errors (blue dashed). The fit is 
at N^LL' order in the R-gap scheme for fti using the central 
values of the correlation matrix given in Eq. ((64}. The center 
of the ellipse are the central values of our final result given in 
Eq. (|68l). 



where the correlation coefficient is significant and reads 

= -0.9176(60) . (65) 

The numbers in the parentheses represent the variance 
from the theory scan. From Eq. (|64p it is straightforward 
to extract the experimental error for as and fii and the 
error due to variations of Q,i and as, respectively: 



exp 



1 



0.0002 , 



cxp 

ill 



CTq^ \Pan\ 



l-pln= 0.009 GcV, 



0.0005, 
0.020 GcV. 



(66) 



For as, the error due to 57i variations is the dominant 
part of the hadronization uncertainty. The blue dashed 
ellipse in Fig. 1151 shows the total error in our final result 
quoted in Eq. (p5)) below. 

The correlation exhibited by the red solid error ellipse 
in Fig.[T5]is indicated by the line describing the semima- 
jor axis 



41.5 GeV 



= 0.1213- a^(mz) ■ 



(67) 



Note that extrapolating this correlation to the extreme 
case where we neglect the nonperturbative corrections 
{^Qi = 0) gives as{mz) — >■ 0.1213. This value is con- 
sistent with the fits in Refs. [l^, [2^ shown in Tab. HI 
which arc dominated by Q = mz where the Monte Carlo 
hadronization uncertainties arc smallest. 
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FIG. 16: Variations of the best fit values for a{rnz) and Sli 
from up (dark shaded blue) and down (light shaded green) 
variations for the theory parameters with respect to the de- 
fault values and in the ranges given in Tab. Illll For the varia- 
tion of the moment Q.2 we use as explained 
in the text. 



Individual Theory Scan Errors 

It is a useful exercise to have a closer look at the size of 
the theory uncertainties caused by the variation of each 
of the theory parameters we vary in our fit procedure 
in order to assess the dominant sources of theory errors. 
In Fig. [TO] two bar charts are shown for the variation of 
the best fit values for as{mz) and ^i{Ra, Pa) at N'^LL' 
order in the R-gap scheme with our default theory pa- 
rameters. The bars show individual up-down variations 
of each of the theory parameters in the ranges given in 
Tab. mil The changes of the best fit values related to 
up variations of the theory parameters are given in dark 
blue and those related to down variations are given in 
light green. 

We see that the dominant theory uncertainties are re- 
lated to variations of the profile functions (ni, ^2, ej, ch) 
and the renormalization scale parameter (n^) for the non- 
singular partonic distribution d(Tns/dr . The uncertain- 
ties related to the numerical errors of the perturbative 
constants (s2, S3, js) as well as the numerical errors in the 
extraction of the nonsingular distribution for small r val- 
ues, (e2, £3) are - with the exception of S2 ~ much smaller 
and do not play an important role. The theory error re- 
lated to the unknown 4-loop contribution to the cusp 
anomalous dimension is negligible. Adding quadratically 
the symmetrized individual errors shown in Fig. [TB] for 
each parameter, we find 0.0006 for as and 0.029 for fii. 
This is about 2/3 of the theoretical uncertainty we have 
obtained by the theory parameter scan, and it demon- 
strates that the theory parameter scan represents a more 
conservative method to estimate the theory error. 

In Fig. [Tni we have also shown the variation of the 
best fit values for as{mz) and ^i{Ra, Pa) due to vari- 
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ations of the second soft function moment parameter 
^2- Our default choice for the parametrization of the 
soft function S™°'^ uses cq = 1 and c„>o ~ with 
A(i?A,/^A) = 0.05 GeV. In this case A is the only vari- 
able parameter of the soft model function 5'™°'^, and fl2 
is predetermined by Eq. (|57p with C2 = 0. As explained 
in Sec. IIVI we modify by setting C2 to nonzero val- 
ues. It is instructive to discuss the Q2 values one should 
consider. From the Cauchy-Schwarz inequality one can 
show that $12/^1 > 1, giving a strict lower bound on 
il2- This bound can only be reached if S™°'^ is a delta- 
function. Moreover, if 5'™°'^ is positive definite, vanishing 
at fc = 0, has a width of order Aqcd, has its maximum at 
a k value of order Aqcd , and has an exponential fall-off 
for large k, then one finds il2/^i < 1.5. We therefore 
adopt the range 1 < il2/^i < 1.5 as a conservative ^2 
variation to carry out an error estimate. For our default 
parametrization we have Q2/^i = 1.18 and changing C2 
between ±0.5 gives a variation of il2/^i between 1.05 
and 1.35. We find that the best fit values for as and fii 
are smooth linear functions of J72/i^i which allows for a 
straightforward extrapolation to the conservative range 
between 1.0 and 1.5. The results for the variations of the 
best fit values for as{mz) and fli for fl2/^i = l.lS^g'^g 

read (<5a,(mz))n, =t"omoll and iSn,)n, =t°oml and 
are also shown in Fig. 1161 The symmetrized version of 
these errors are included in our final results. For our final 
results for as{mz) we add the uncertainties from J7i and 
the one from f22 quadratically giving the total hadroniza- 
tion error. For f^i(i?A, /^a) we quote the error due to f22 
separately. 

Final Results 

As our final result for as{mz) and f2i(i?A, A^a), obtained 
at N'^LL' order in the R-gap scheme for fii, including 
bottom quark mass and QED corrections we obtain 

as{mz) = 0.1135 ± (0.0002)exp 

± (0.0005)hadr ± (0.0009)pcrt, 

l^i(i?A,MA) = 0.323 ± (0.009)cxp ± (0.013)a, 

± (0.020)„^(„,,) ± (0.045)pert GcV, (68) 

where i?A = = 2 GeV and we quote individual 1- 
sigma errors for each parameter. Eq. (j68p is the main 
result of this work. In Fig. [15] (blue dashed line) and 
Fig. [TTk (thick dark red line) we have displayed the cor- 
responding combined total (experimental-|-theoretical) 
standard error ellipse. To obtain the combined ellipse we 
take the theory uncertainties given in Tabs. |IV] and |V] to- 
gether with the f22 uncertainties, adding them in quadra- 
ture. The central values in Eq. are determined by 
the average of the respective maximal and minimal val- 
ues of the theory scan, and are very close to the central 
values obtained when running with our default theory 
parameters. The fit has /dof ~ 0.91 with a variation 
of ±0.03 for the displayed scan points. Having added the 



theory scan and f22 uncertainties reduces the correlation 
coefficient in Eq. ^ to p^°}f = -0.212. As a compar- 
ison we have also shown in Fig. [TTb the combined total 
(experimental+theoretical) error ellipse at N^LL' in the 
MS scheme for f2i where the ©(Aqcd) renormalon is not 
subtracted. 

Since our treatment of the correlation of the system- 
atic experimental errors is based on the minimal over- 
lap model, it is instructive to also examine the results 
treating all the systematic experimental errors as uncor- 
rclated. At N^LL' order in the R-gap scheme the re- 
sults that arc analogous to Eqs. (|68|) read as{mz) = 
0.1141 ± (0.0002)cxp ± (0.0005)hadr ± (O.OOlO)pcrt and 
f^i(i?A,/iA) = 0.303±(0.006)exp±(0.013)n,±(0.022)„^± 
(0.055)pcrt GeV with a combined correlation coefficient of 
PqO^' ~ —0.180. The results arc compatible with the re- 
sults of Eqs. and indicate that the ignorance of the 
exact correlation of the systematic experimental errors 
does not crucially affect the outcome of the fit. 

Data Set Choice 

We now address the question to which extent the results 
of Eqs. depend on the thrust ranges contained in the 
global data set used for the fits. Our default global data 
set accounts for all experimental thrust bins for Q > 35 
in the intervals [TminjTmax] = [6/(5,0.33]. (See Sec. IVll 
for more details.) This default global data set is the 
outcome of a compromise that (i) keeps the r interval 
large to increase statistics, (ii) sets Tmin sufficiently large 
such that the impact of the soft function moments 17^ 
with i > 2 is small and (iii) takes Tmax sufficiently low 
to exclude the far-tail region where the missing order 
Q^s Aqcd IQ corrections potentially become important. 

In Fig. [T7] the best fits and the respective experimen- 
tal 39% and 68% CL error ellipses for the default values 
of the theory parameters given in Tab. Illll are shown for 
global data sets based on different t intervals. The re- 
sults for the various t intervals are each given in different 
colors. The results for our default global data set is given 
in red color, and the subscript "strict" for some intervals 
means that bins are included in the data set if more than 
half their range is contained within the interval. For in- 
tervals without a subscript the criterion for selecting bins 
close to the boundaries of the r interval is less strict and 
gencrically, if the Tmin and Tmax values fall in such bins, 
these bins are included. The numbers in superscript for 
each of the r intervals given in the figure refers to the to- 
tal number of bins contained in the global data set. We 
observe that the main effect on the outcome of the fit 
is related to the choice of r„iin and to the total number 
of bins. Interestingly all error ellipses have very similar 
correlation and are lined up approximately along the line 



50.2 GeV 



= 0.1200- as(mz) 



(69) 



Lowering T^i^ increases the dependence on ^2 and leads 
to smaller as and larger fii values. On the other hand, 
increasing Tmin leads to a smaller data set and to larger 
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as{mz) 

FIG. 17: The smaller elongated ellipses show the experimental 
39% CL error (1-sigma for Os) and best fit points for different 
global data sets at N"^LL' order in the R-gap scheme and 
including bottom quark mass and QED effects. The default 
theory parameters given in Tab. lIIII are employed. The larger 
ellipses show the combined theoretical plus experimental error 
for our default data set with 39% CL (solid, 1-sigma for one 
dimension) and 68% CL (dashed). 



experimental error ellipses, henee to larger uneertainties. 

It is an interesting but expeeted outcome of the fits 
that the pure experimental error for as (the uncertainty 
of as for fixed central fii) depends fairly weakly on the 
T range and the size of the global data sets shown in 
Fig. [iTl If we had a perfect theory description then we 
would expect that the centers and the sizes of the error 
ellipses would be statistically compatible. Here this is 
not the case, and one should interpret the spread of the 
ellipses shown in Fig. [17] as being related to the theo- 
retical uncertainty contained in our N'^LL' order predic- 
tions. In Fig. [17] we have also displayed the combined 
(experimental and theoretical) 39% CL standard error 
ellipse from our default global data set which was al- 
ready shown in Fig. [TTb (and is 1-sigma, 68% CL, for 
either one dimensional projection). We also show the 
68% CL error ellipse by a dashed red line, which corre- 
sponds to 1-sigma knowledge for both parameters. As 
we have shown above, the error in both the dashed and 
solid larger ellipses is dominated by the theory scan un- 
certainties, see Eqs. (|68|) . The spread of the error ellipses 
from the different global data sets is compatible with the 
1-sigma interpretation of our theoretical error estimate, 
and hence is already represented in our final results. 

Analysis without Power Corrections 

Using the simple assumption that the thrust distribution 
in the tail region is proportional to as and that the main 





Q!s(?nz)±(pert. error) 


xV(dof) 


N^LL' with 


0.1135 ±0.0009 


0.91 


N^LL' with Clf^ 


0.1146 ±0.0021 


1.00 


N^LL' without Sr""^ 


0.1241 ±0.0034 


1.26 


0{al) fixed-order 
without 5'"°'* 


0.1295 ±0.0046 


1.12 



TABLE VII: Comparison of global fit results for our full anal- 
ysis to a fit where the renormalon is not canceled with Qi, a 
fit without 5™°'* (meaning without power corrections with 
S™°'^(fc) — S{k)), and a fit at fixed order without power cor- 
rections and log resummation. All results include bottom 
mass and QED corrections. 



effect of power corrections is a shift of the distribution 
in T, we have estimated in Sec. [T] that a 300 MeV power 
correction will lead to an extraction of as from Q = niz 
data that is Sag /as — (—9 ± 3)% lower than an anal- 
ysis without power corrections. In our theory code we 
can easily eliminate all nonperturbative effects by set- 
ting S'^°'^{k) = 5{k) and A = 5 = 0. At N^LL' or- 
der and using our scan method to determine the per- 
turbative uncertainty a global fit to our default data set 
yields as{mz) = 0.1241 ± (0.0034)port which is indeed 
9% larger than our main result in Eq. (|68p which ac- 
counts for nonperturbative effects. It is also interesting 
to do the same fit with a purely fixed-order code, which 
we can do by setting ns = ^^J = fJ-H to eliminate the 
summation of logarithms. The corresponding fit yields 
C(s{iTiz) = 0.1295±(0.0046)port, where the displayed error 
has again been determined from the theory scan which in 
this case accounts for variations of fiH and the numerical 
uncertainties associated with £2 and £3. (A comparison 
with Ref. [l^] is given below in Sec. IIXI ) 

These results have been collected in Tab. IVIll together 
with the as results of our analyses with power corrections 
in the R-gap and the MS schemes. For completeness we 
have also displayed the respective x^/dof values which 
were determined by the average of the maximal and the 
minimum values obtained in the scan. 



VIII. FAR- TAIL AND PEAK PREDICTIONS 

The factorization formula ([4]) can be simultaneously used 
in the peak, tail, and far-tail regions. To conclude the 
discussion of the numerical results of our global analysis 
in the tail region, we use the results obtained from this 
tail fit to make predictions in the peak and the far-tail 
regions. 

In Fig. [18] we compare predictions from our full N'^LL' 
code in the R-gap scheme (solid red line) to the accurate 
ALEPH data aX Q ~ mz in the far-tail region. As input 
for as{mz) and fli we use our main result of Eq. (|68)) 
and all other theory parameters are set to their default 
values (sec Tab. lIII|) . We find excellent agreement within 
the theoretical uncertainties (pink band). Key features 
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FIG. 18: Thrust distributions in the far-tail region at N^LL' 
order with QED and mi, corrections included at Q = mz to- 
gether with data from ALEPH. The red solid line is the cross 
section in the R-gap scheme using asijnz) and fii obtained 
from fits using our full code, see Eq. (|68p . The hght red band 
is the perturbative uncertainty obtained from the theory scan 
method. The red dashed line shows the distribution with the 
same but without power corrections. The Ught solid blue 
line shows the result of a full N'^LL' fit with the BS profile 
that does not properly treat the multijet thresholds. The 
short dashed green line shows predictions at N^LL' with the 
BS profile, without power corrections, and with the value of 
as{mz) obtained from the fit in Ref. [23|. All theory results 
are binned in the same manner as the experimental data, and 
then connected by lines. 



of our theoretical result in Eq. (j?]) that arc important in 
this far-tail region are i) the nonperturbative correction 
from r^i, and ii) the merging of fisiT), /ij(r), and fin 
toward fis = /^j = [i-H at r = 0.5 in the profile func- 
tions, which properly treats the cancellations occurring 
at multijet thresholds. To illustrate the importance of 
ill we show the long-dashed red line in Fig. [15] which has 
the same value of cxs(mz)., but turns off the nonpertur- 
bative corrections. To illustrate the importance of the 
treatment of multijet thresholds in our profile function, 
we take the BS profile which does not account for the 
thresholds (the BS profile is defined and discussed below 
in Sec. HX)) . and use the smaller asimz) and larger fii 
that are obtained from the global fit in this case. The 
result is shown by the solid light blue line in Fig. [T51 
which begins to deviate from the data for r > 0.36 and 
gives a cross section that docs not fall to zero at r = 0.5. 
The fact that as{mz) is smaller by 0.0034 for the light 
blue line, relative to the solid red line, indicates that the 
proper theoretical description of the cross section in the 
far-tail region has an important impact on the fit done 
in the tail region. The final curve shown in Fig. [18] is the 
short-dashed green line, which is the result at the level 
of precision of the analysis by Becher and Schwartz in 
Ref. It uses the BS profile, has no power corrcc- 



FIG. 19: Thrust cross section for the result of the N^LL' fit, 
with QED and rUb corrections included at Q = mz- The 
red solid line is the cross section in the R-gap scheme using 
as{mz) and Sli obtained from fits using our full code, see 
Eq. (|68p . The red dashed line shows the distribution with the 
same as but without power corrections. The short-dashed 
green line shows predictions at N^LL' with the BS profile, 
without power corrections, and with the value of as{mz) ob- 
tained from the fit in Ref. Data from ALEPH, DELPHI, 
L3, SLD, and OPAL are also shown. 



tions, and has the value of obtained from the fit in 
Ref. [1^. It also misses the Q — mz data in this re- 
gion. The results of other 0{al) thrust analyses, such as 
Davison and Webber and Dissertori et al. [22l. [25j. 
significantly undershoot the data in this far-tail region. 
To the best of our knowledge, the theoretical cross sec- 
tion presented here is the first to obtain predictions in 
this far-tail region that agree with the data. Note that 
our analysis does include some ©(a^ Aqcd/Q) power cor- 
rections through the use of Eq. ([^^ . It does not account 
for the full set of 0(q;sAqcd/Q) power corrections as 
indicated in Eq. ([4]) (see also Tab. Illb). but the agree- 
ment with the experimental data seems to indicate that 
missing power corrections may be smaller than expected. 

Unbinned predictions for the thrust cross section at 
Q = mz in the peak region arc shown in Fig. 1191 The 
green dashed curve shows the result at the level of pre- 
cision in Becher and Schwartz, that is N'^LL', with the 
BS profile, without power corrections, and with the value 
of as{mz) = 0.1172 obtained from their fit. This purely 
perturbative result peaks to the left of the data. With 
the smaller value of asimz) obtained from our fit, the 
result with no power corrections peaks even slightly fur- 
ther to the left, as shown by the long-dashed red curve. 
In contrast, the red solid curve shows the prediction from 



See the top panel of Fig. 9 in Ref. [531, the top left panel of Fig. 4 
in Ref. [23], and the left panel of Fig. 2 in Ref. [2511 . 



our full N^LL' code in the R-gap scheme with our central 
fit values of as{mz) and fli given in Eq. It clearly 

indicates that the value of fli obtained from the fit in 
the tail region shifts the theory prediction in the peak 
region much closer to the experimental data. The resid- 
ual difference between the solid red theory curve and the 
experimental data can be attributed to the fact that the 
peak is sensitive to power corrections from higher mo- 
ments, rifc>2, which have not been fit in our analysis. In 
our theoretical cross section result this would correspond 
to fitting A(i?A,A*A), and a subset of the higher coeffi- 
cients Ci>i. The Ci>i were all set to zero in the curves 
shown here, and we leave the presentation of results of 
this extended fit to a future publication. 



IX. CROSS CHECKS AND COMPARISONS 

The result for as{mz) we obtain from our global 
N-^LL' analysis in the R-gap scheme with 487 bins 
given in Eq. is consistent at 1-sigma with the re- 
suh of Davison and Webber [H] (as(TOz) = 0.1164 ± 
(0.0022)hadr-Hexp ± (0.0017)port)- They also carried out 
a global thrust analysis with a total of 430 experimen- 
tal bins. As explained in Sec. HI in their theory formula 
nonperturbative effects were included as a power correc- 
tion in the effective coupling model which was fit from 
the experimental data, and their approach also accounts 
for a renormalon subtraction of the perturbative distri- 
bution. In these respects their analysis is similar to ours. 
However, it differs as their theory formula contains only 
resummation of logarithms at NLL order, and it also 
uses a different renormalon subtraction scheme which is 
based on the running coupling approximation for the sub- 
traction corrections and does not account for the resum- 
mation of large logarithms. Moreover the separation of 
singular and nonsingular perturbative contributions and 
method to turn off the log resummation at large t is 
not equivalent to the one we employ. The difference be- 
tween their central value and perturbative error and our 
Eq. ([551) can be attributed to these items. Their com- 
bined hadronization and experimental uncertainty uti- 
lizes an error rescahng using the value Xmin/'^'^^ ~ 1-09 
obtained for their best fit. 

On the other hand, our main result for as{mz) given in 
Eq. is smaller than the results of Dissertori et al. [l^l 
by 2.9-sigma, of Dissertori et al. [25| by 2.2-sigma, and 
of Becher and Schwartz [2^ by 1.6-sigma. (These results 
are displayed in Table HI) In these analyses as{mz) was 
determined from fits to data for individual Q values and, 
as explained in Sec. [H nonperturbative corrections and 
their associated uncertainty were taken from Monte Carlo 
generators in Dissertori et al., or left out from the fit 
and used to assign the hadronization uncertainty for the 
final result in Becher and Schwartz. It is possible to 
turn off pieces of our theoretical code to reproduce the 
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j_jauc;i iiiiciiij 


rjiiLi gy 


BS 
results [20] 


our BS 
profile 


default 
profile 


ALEPH 


91.2 GeV 


0.1168(1) 


0.1170 


0.1223 


ALEPH 


133 GeV 


0.1183(37) 


0.1187 


0.1235 


ALEPH 


161 GeV 


0.1263(70) 


0.1270 


0.1328 


ALEPH 


172 GeV 


0.1059(80) 


0.1060 


0.1088 


ALEPH 


183 GeV 


0.1160(43) 


0.1166 


0.1205 


ALEPH 


189 GeV 


0.1203(22) 


0.1214 


0.1260 


ALEPH 


200 GeV 


0.1175(23) 


0.1182 


0.1224 


ALEPH 


206 GeV 


0.1140(23) 


0.1149 


0.1185 


OPAL 


91 GeV 


0.1189(1) 


0.1198 


0.1251 


OPAL 


133 GeV 


0.1165(38) 


0.1175 


0.1218 


OPAL 


177 GeV 


0.1153(33) 


0.1160 


0.1200 


OPAL 


197 GeV 


0.1189(14) 


0.1197 


0.1241 


average 




0.1172(10) 


0.1180 


0.1221 


global tit 

(stat) 
global fit 
(stat+syst) 


all Q 
all Q 




0.1188 
0.1192 


0.1242 
0.1245 



TABLE VIII: Comparison of the results for Qs(mz) quoted by 
Becher and Schwartz in Ref. [2^ with results we obtain from 
our adapted code where power corrections, the mj, and QED 
corrections, the 0{a'i) axial singlet corrections are neglected. 
The 0(al) nonlogarithmic constants hs and S3 are set to the 
values used in Ref. [2011 as described in the text. We follow the 
fit approach of Ref. [20| and employ their profile functions for 
the nonsingular, hard, jet and soft scales, with results shown 
in the column labeled "our BS profile" . In the last column we 
show results with this same code, but using our default profile 
functions. The errors in the third column are the statistical 
experimental uncertainty. 



perturbative precision of the codes used in Refs. [22j^^ 
and poj . It is the main purpose of the remainder of this 
section to show the outcome of the fits based on these 
modified theory codes. We show in particular, that the 
main reason why the above results for asijnz) arc higher 
than our result of Eq. ([S5[) is related to the fact that the 
nonperturbative corrections extracted from Monte Carlo 
generators at Q = mz are substantially smaller than and 
incompatible with the ones obtained from our fit of the 
field theory power correction parameter fii. In SecUwe 
already discussed why the use of Monte Carlo generators 
to estimate nonperturbative corrections in high-precision 
perturbative predictions is problematic. 

We start with an examination related to the code used 
by Becher and Schwartz which has N'^LL' accuracy 
but docs not include power corrections or renormalon 
subtractions. This treatment can be reproduced in our 
factorization formula by turning off the nonperturbative 



We do not attempt to reproduce the NLL/Ofg^) code of Ref. [25| 
as the final outcome is similar to Ref. |22|| . 
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soft nonperturbative function by setting S™°'^{k) = 6{k) 
and A = 6 = 0. Moreover they used the central scale 
setting fiH ~ Q, — QVt' and /ig = Qt. We can repro- 
duce this from our profile functions for fiQ = ni = ej = 0, 
t2 = 3/2 and en ^ ng = 1, which we call the BS profile 
setting. The BS profile functions for hj{t) and fisiT) 
are shown by dashed curves in Fig. [5] (Note that the 
BS profile setting docs not cause fis-, fJ-J^ a-nd fin to 
merge in the far-tail region and become equal at r — 0.5, 
which is needed to switch off the SCET resummation 
of logarithms in the multijet region to satisfy the con- 
straints from multijet thresholds.) Becher and Schwartz 
set the 0{al) nonlogarithmic correction in the Euclidean 
hard factor C{—q^) to zero (with Hq ~ 1^(5^)1^ for 
= Q'^ > 0), which in our notation corresponds to 
/i3 ~ 11771.50 (somewhat larger than the now known 
hs). We also set sa = -40.1 (see Ref. [13 H^) and 
S3 = -324.631 for the non-logarithmic O(a^) and 0{al) 
constants in the soft function (both within our range of 
uncertainties). The value for S3 corresponds to setting 
the 0{a^) nonlogarithmic corrections in the expanded 
position space soft function to zero. Finally, we also turn 
off our QED and bottom quark mass corrections and the 
0{aj.) axial singlet corrections, and use the fixed-order 
normalization from Eq. ([5^ . For the fit procedure we 
follow Becher and Schwartz and analyze all ALEPH and 
OPAL data for individual Q values in the r ranges given 
in their work and account only for statistical experimen- 
tal errors in the functions. The outcome of the fits for 
as{mz) at N^LL' order is given in the fourth column of 
Tab. I Vim The third column shows their central values 
and the respective statistical experimental errors as given 
in Ref. The numbers we obtain are 0.0001 to 0.0011 
higher than their central values, and we attribute this 
discrepancy to the nonsingular contributions.^^ (Becher 
and Schwartz also used a difference of cumulants for their 
fits, as in Eq. (|62p with the choice fi = ti and fa = T2, 
rather than integrating dcr/dr as we do for the table. The 
spurious contribution induced by this choice has a signifi- 
cant effect on the values, but a small effect on as{mz), 
changing the values shown in the table by < 0.0003. For 
cumulants that use fi = fa = (ti+T2)/2 with no spurious 
contribution, the difference from our integrated distribu- 
tion results is reduced to < 0.0001 for as{mz), and 
values arc almost unaffected.) 

The numbers obtained at N'^LL' above are significantly 
larger than our central fit result as{mz) — 0.1135 shown 
in Eq. obtained from our full code. These differ- 



ences are mainly related to the nonperturbative power 
correction and partly due to the BS profile setting. To 
distinguish these two and other effects we can take the 
purely perturbative code described above and turn back 
to our default setting for the profile functions with the 
parameters given in Tab. IIIII The results are shown in 
the fifth column in Tab. IVIIII using again only statisti- 
cal experimental errors in the x^ functions. The (Xsimz) 
values using our default profile functions are by 0.0028 
to 0.0058 larger than for the BS profile setting in the 
fourth column. (The fifth column results again inte- 
grate the distribution over each bin rather than using 
differences of cumulants, which for our profile is impor- 
tant for the reasons discussed in Sec. |Vl^^) A similar 
difference arises from a global fit to our default data set 
of Sec. IVII using the same fit procedure: For the BS pro- 
file setting we obtain as{mz) = 0.1189, while the default 
profile setting gives as{mz) = 0.1242 (second to last line 
of Tab. IVIIip . Using instead the x^-analysis of our main 
analysis which includes the experimental systematical er- 
rors we obtain as{mz) = 0.1192 for the BS profile setting 
and as{mz) = 0.1245 for the default profile setting (last 
line of Tab. NTU^ . The latter resuft is by 0.0110 larger 
than our 0.1135 central fit resuft in Eq. §8^. This 10% 
effect is almost entirely coming from the power correc- 
tion r^i. The difference of 0.3% to the full perturbative 
result of as{mz) = 0.1241 given in Table IVlIl illustrates 
the combined effect of the QED, the bottom quark mass 
and the 0(0^) axial singlet corrections and the 0{al) 
hard constant h^. 

Finally, let us examine the results related to the code 
used by Dissertori et al. in Ref. [13], which uses the 
fixed-order 0{a^) results without a resummation of loga- 
rithms, but accounts for nonperturbative corrections de- 
termined from the difference of running Monte Carlo gen- 
erators in parton and hadron level modes. Since in this 
work we are not concerned with extracting the parton- 
hadron level transfer matrix from Monte Carlo genera- 
tors, we use in the following our code neglecting power 
corrections by setting S^°^{k) = S{k), setting A = 6 = 0, 
and setting hh = y^j = pi-s- The latter switches off the 
log resummation factors in Eq. (jl]) such that only the 
0{a'l) fixed order expression remains. We also include 
the mb corrections, but neglect QED effects. Since these 
modifications give us a code that does not contain non- 
perturbative corrections, the differences to Ref. [13] we 
obtain will serve as a quantitative illustration for the size 



^ Becher and Schwartz uncovered a numerical problem with the 
original EERAD3 code at very small t, which correspondingly 
had an impact on the nonsingular function used in their analysis 
which was extracted from EERAD3. When their nonsingular 
distribution is updated to results from the new EERAD3 code 
they become significantly closer to ours, differing by < 0.0002. 
We thank M. Schwartz for correspondence about this and for 
providing us with their new fit values. 



With our full code, which accounts in particular for power correc- 
tions and rcnormalon subtractions, the shift due to the modified 
profile functions becomes smaller; shifts in as(mz) of 0.005 be- 
come 0.003. 

Using the cumulant method with ti = ri and f 2 = T2 in Eq. II62I I , 
which has a spurious contribution, changes the values in the fifth 
column of Tab. |llT] by about -0.003 to -0.005. On the other 
hand, using the cumulant method without a spurious contribu- 
tion, ri = f 2 = (ti -|-r2)/2, changes the values in the fifth column 
by < 0.0001. 
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Dissertori et al. 
results [22] 


Our fixed 
order code 


ALEPH 


9L2GeV 


0.1274(3) 


0.1281 


ALEPH 


133 GeV 


0.1197(35) 


0.1289 


ALEPH 


161 GeV 


0.1239(54) 


0.1391 


ALEPH 


172 GeV 


0.1101(72) 


0.1117 




ioo Oe V 


U.iioz^oz ) 


U.iz4( 


ALEPH 


189 GeV 


0.1140(20) 


0.1295 


ALEPH 


200 GeV 


0.1094(22) 


0.1260 


ALEPH 


206 GeV 


0.1075(21) 


0.1214 



TABLE IX: Comparison of the thrust results quoted in 
Ref. [l^l with our numerical reproduction. For this numerical 
exercise we have used their procedure to get the error matrix 
for the experimental data. This amounts to considering only 
the statistical errors in an uncorrelated way, with the resulting 
experimental error shown in the third column. Whereas in the 
code of Ref. [S^l hadronization corrections are included deter- 
mined from Monte Carlo simulations our numbers are based 
on a pure partonic code neglecting nonperturbative effects. 
We use the default value for the scale setting, i.e. jj, — Q. 



of the hadronization corrections obtained by a transfer 
matrix from the Monte Carlo generators PYTHIA, HER- 
WIG, and ARIADNE, tuned to global hadronic obscrv- 
ables at mz- 



For the fits for as{mz) we follow Dissertori et al. [22 1 
analyzing ALEPH data for individual Q values in the r 
ranges given in their work and accounting only for sta- 
tistical experimental errors in the functions. The re- 
sults of Dissertori et al. and the outcome for our best 
fits are given in the third and fourth column of Tab. IIXI 
respectively. We have also quoted the respective statis- 
tical errors from Ref. [l^. For the high statistics data 
at Q = mz our as{mz) result is larger than theirs, 
but the discrepancy amounts to only 0.0007 which is a 
0.5% shift in as{mz). This illustrates the small size of 
the nonperturbative hadronization corrections encoded 
in the Monte Carlo transfer matrix at Q = niz- This 
is clearly incompatible with the size of the nonperturba- 
tive correction we have obtained from simultaneous fits 
of as and fii, confirming the concerns on Monte Carlo 
hadronization corrections explained in Sec. HI Interest- 
ingly, with the exception of Q = 172 GeV, our fixed- 
order results for all Q are relatively stable and close to 
the result at Q = mz, while their as{mz) values, which 
use the transfer matrix for nonperturbative effects, are 
systematically lower for Q > mz by 7 to 13%. Thus 
the nonperturbative effects from the Monte Carlo trans- 
fer matrix are substantially larger for Q > mz-'^^ The 



same behavior is also visible in the results of Ref. [25|, 
which includes NLL resummation of logarithms. Since 
the transfer matrix is obtained from Monte Carlo tuned 
to the more accurate Q = mz data, we believe that 
this issue deserves further investigation. To complete 
the discussion we use the same fixed-order theory code 
to quote results for a global fit to our default data set. 
Using the fit procedure as described in Sec. IVII we obtain 
as{mz) = 0.1300 ± (0.0047)port- (The corresponding er- 
rors obtained from the error band method are given in 
the fourth Une of Tab. Ell) 



X. CONCLUSIONS 

In this work we have provided a factorization formula 
for the thrust distribution in e^e~ annihilation which in- 
corporates the previously known 0{al) and ©(a^) per- 
turbative QCD corrections and summation of large loga- 
rithms at N'^LL order for the singular terms in the dijet 
limit where the thrust variable r = 1 — T is small. The 
factorization formula used here incorporates a systematic 
description of nonperturbative effects with a soft function 
defined in field theory. The soft function describes the 
dynamics of soft particle radiation at large angles. We 
have also accounted for bottom mass and QED photon 
effects for fixed-order contributions as well as for the sum- 
mation of QED logarithms. With specifically designed r- 
dependent profile functions for the renormalization scales 
the factorization formula can be applied in the peak, tail 
and far-tail regions of the thrust distribution. It has 
all nonperturbative effects accounted for up to terms of 
©(asAqcD/Q), which is parametrically smaller than the 
remaining perturbative uncertainty (< 2% for Q = mz) 
of the thrust distribution predictions in the tail region 
where we carried out the fits to the experimental data. 

In the tail region, 2Aqcd/Q ^ 1/3, the domi- 

nant effects of the nonperturbative soft function are en- 
coded in its first moment f2i, which is a power correc- 
tion to the cross section. Fitting to tail data at multiple 
Qs as we did in this work, the strong coupling as{mz) 
and the moment ili can be simultaneously determined. 
An essential ingredient to reduce the theoretical uncer- 
tainties to the level of < 2% in the thrust distribution 
is our use of a short-distance scheme for Oi, called the 
R-gap scheme, that induces subtractions related to an 
C(Aqcd) renormalon contained in the MS perturbative 
thrust cross section from large angle soft gluon radia- 
tion. The R-gap scheme introduces an additional scale 
that leads to large logarithms in the subtractions, and 
we carry out a summation of these additional logarithms 
with renormalization group equations in the variable R. 
The R-gap scheme reduces the perturbative uncertainties 
in our best highest order theory code by roughly a factor 



■^^ Note that the weighted average of the Q > mz thrust results of 
Dissertori et al. is as{mz) = 0.1121 and is consistent with our 
result in Eq. I|68| l within the larger uncertainties. Also note that 
the Q dependence of our Qi{R, R)/Q power correction is affected 



by its anomalous dimension, cf. Fig. [6] 
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FIG. 20: Comparison of selected determinations of as{mz) defined in the MS scheme. 



of two compared to the pure MS definition, f2i, where 
renormalon effects are not treated. 

The code we use in this analysis represents the most 
complete theoretical treatment of thrust existing at this 
time. As our final result we obtain 



as{mz) = 0.1135 ± 0.0011, 
f^i(i?A,MA) = 0.323 ± 0.051 GcV, 



(70) 



where as is defined in the MS scheme, and ili in the R- 
gap scheme at the reference scales i?A = = 2 GcV. 
Here the respective total 1-sigma errors are shown. The 
results with individual 1-sigma errors quoted separately 
for the different sources of uncertainties are given in 
Eq. Neglecting the nonperturbative effects incor- 

porated in the soft function, and in particular from 
the fits gives as{mz) = 0.1241 which exceeds the result 
in Eq. ([70)) by 9%. This is consistent with a simple scal- 
ing argument one can derive from experimental data, see 
Eq. dS]) in Sec. H 

Analyses of event shapes with a simultaneous fit of 
as and a power correction have been carried out earlier 
with the effective coupling model. Davison and Web- 
ber [23| analyzed the thrust distribution and determined 
asimz) = 0.1164 ± 0.0028 also using 0{al) fixed-order 
input, but implementing the summation of logarithms 
only at NLL order (for further discussion see Sec. lIXp . 
Recently Gehrmann et al. [95| analyzed moments of dif- 
ferent event shape distributions, also with the effective 
coupling model, and obtained as{mz) = 0.1153 ± 0.0029 
using fixed-order perturbation theory at O(a^). Both 
analyses neglected bottom mass and QED corrections. 
Our result in Eq. ([70)) is compatible with these analyses 
at 1-sigma, but has smaller uncertainties. 

These results and our result for as{mz) in Eq. ()70p 
are substantially smaller than the results of event shape 
analyses employing input from Monte Carlo generators 



to determine nonperturbative effects. We emphasize that 
using parton-to-hadron level transfer matrices obtained 
from Monte Carlo generators to incorporate nonpertur- 
bative effects is not compatible with a high-order theo- 
retical analysis such as ours, and thus analyses relying on 
such Monte Carlo input contain systematic errors in the 
determination of as from thrust data. The small effect 
of hadronization corrections on thrust observed in Monte 
Carlo generators at Q = mz and the corresponding small 
shift in asimz) do not agree with the 9% shift we have 
obtained from our fits as mentioned above. For the rea- 
sons discussed earlier, we believe Monte Carlo should not 
be used for hadronization uncertainties in higher order 
analyses. 

Although our theoretical approach represents the most 
complete treatment of thrust at this time, and all sources 
of uncertainties known to us have been incorporated in 
our error budget, there are a number of theoretical is- 
sues related to subleading contributions that deserve fur- 
ther investigation. These issues include (i) the summa- 
tion of logarithms for the nonsingular partonic cross sec- 
tion, (ii) the structure of the ©(asAqco/Q) power cor- 
rections, (iii) analytic perturbative computations of the 
0{a^) and O(a^) nonlogarithmic coefficients S2 and S3 
in the partonic soft function, the Oia^) nonlogarithmic 
coefficient ja in the partonic jet function, and the 4-loop 
QCD cusp anomalous dimension Pg"''''. Concerning is- 
sue (i) we have incorporated in our analysis the non- 
singular contributions in fixed-order perturbation theory 
and estimated the uncertainty related to the higher order 
logarithms through the usual rcnormalization scale vari- 
ation. Further theoretical work is needed to derive the 
rcnormalization group structure of subleading jet, soft, 
and hard functions in the nonsingular contributions and 
to use these results to sum the corresponding logarithms. 
Concerning issue (ii) we have shown that our theoretical 
description for the thrust distribution contains a remain- 
ing theoretical uncertainty from nonperturbative effects 
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of order ©(agAqcD/Q)- Parametrically, this uncertainty 
is substantially smaller than the perturbative error of 
about 1.7% for the thrust distribution in the tail region 
at LEP-I energies that is contained in our best theory 
code. Furthermore, our predictions in the far-tail region 
at Q = mz appear to indicate that the dominant cor- 
rections of this order are already captured in our setup. 
Nevertheless a systematic analysis of these subleading ef- 
fects is certainly warranted. 

Apart from investigating these theoretical issues, it is 
also warranted to apply the high-precision approach us- 
ing soft-collinear effective theory to other event shape 
distributions in order to validate the result in Eq. ([70| . 
Event shapes that can be clearly treated with simi- 
lar techniques are: heavy and li ght jet masses, the C- 
parameter, and the angularities [96l |97|. For many of 
these event shapes it has been proven field theoretically 
that the same parameter fii describes the leading power 
corrections in the tail region , although there might be 
caveats related to the experimental treatment of hadron 
masses [13, [sl]. Thus, one has the potential to extend the 
analysis done here to include additional data without ad- 
ditional parameters. An analysis for the heavy jet mass 
accounting for perturbative contributions at N'^LL in MS 
with different profile functions and a simple soft function 
model for power corrections without renormalon subtrac- 
tions, was recently carried out in Ref. j99| . providing a 
first step in this direction. 

To conclude this work we cannot resist comparing our 
result for as (jnz) with the results of a selection of anal- 
yses using other techniques and observables, as shown 
in Fig. [201 We include a N'^LO analysis of da ta fr om 
deep inelastic scattering in the nonsinglet channel jlOOf ^^. 
the recent HPQCD lattice determination based o n fit - 
ting Wilson loops and the T-T' mass difference |l05l |. 



the result from fits to electroweak precision observables 
based on the Gfitter p ackag e |l06| . analyses of r-decay 
data using fixed-o rder [107| and contour-improved per- 
turbation theory Il08| . together with an average of r 
results from Ref. |l09| . Final ly we als o show a collec- 
tion of Qfs-averages from Refs. 109l4lll |. The DIS result 
is consistent with our fit result, whereas the deviation 
from HPQCD is 3.5cr. It is interesting to note that the 
high energy extractions from thrust and DIS appear to 
be smaller than the low energy extractions from Lattice 
and T decays. 
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Appendix A: Formulae 

In this appendix we collect all the remaining formulas 
used in our analysis for the case of massless quarks. The 
total hadronic cross section at tree level at the energies 
we are considering is 

MQ)^ E [^LW) + <c(Q)], (Al) 

where Q is the cm. energy. For a quark of flavor q the 
tree level axial-vector and vector cross sections are 



qa 



- Nr. 



Q\vl + al)al 



Q2)2 



^ z 



(A2) 



qv 



47ra^ 



2 eqVqVe ' 



(m| - Q2)2 + 



^ z 



nvi+aiH 



(m 



where eq is the electric charge of the quark, and 



Tl -2eq sin^ 



rpq 
-^3 



sin(2 9iv) 



sin(2 9w) 



(A3) 



Analyses studying as with data tha t depends also on the gluon 
PDF have been carried out in Refs. [lO]Hl04ll . 



Here Tg is the third component of the weak isospin, and 
9w is the weak mixing angle. For our numerics we use 
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the following values: 

sin^ 6'vi/ = 0.23119 , = 91.187 GeV , 

Tz = 2.4952 GeV, TOt = 172GeV, 
mfc = 4.2GeV, a(mz) 1/127.925 . (A4) 

Singular Cross Section Formula 

To simplify the numerical evaluation of the singular part 
of the differential cross section given in Eq. pT|) we take 
fi = so that U}{s — s' , jjLj, /ij) = 5{s — s') and express 
the result in the following form 



j dkP{Q,QT-k,^ij) 

e-2«(^^''^)^ ^mod A _ 2A(i?, iis)) , (A5) 



where the perturbative corrections from the par- 
tonic soft function, jet function, and soft evo- 
lution factor are contained in P{Q^ fc, A' j) ~ 

/ds/ dk' Jr{s,^lJ)Ul{k',^iJ,^ls)s^'''\k-k' -s/Q,^is)■ 

The integrals in P can be carried out explicitly so 
that it is given by a simple set of functions. The soft 
nonperturbative function Sf°^{k — 2 A) is discussed in 
Sec. HVl and in Eq. (|A5|) we have integrated by parts 
so the derivative in the exponential with the S{R,fis) 
acts on this nonperturbative function. Hq, Jr, S'p^''* 
and exp{—2S{R, fj,s)d/dk) (cf. Eq. ^7^) involve series 
in as(/i/j), Q!s(a*j), and as{ns) with no large logs, and 
in our numerical analysis we expand the product of 
these series out, order- by-order in as- This expansion is 
crucial for S^'^'^^{k, fis) and cxp{—26{R, fis)d/dk) since 
it is needed to allow the rcnormalon in the two series to 
cancel. 

For simplicity where possible we give ingredients in a 
numerical form for SU(3) color with Uf = 5 active flavors. 
The vector hard function to 0{al) is [H,!!!-!!! 



HQiQ,fi) 

= 1 + a^(Ath) (0.745808- 1.27324Lq-0.848826L^ 
-t- ^2.27587 - 0.0251035 Lq - 1.06592 

+ 0.735517L| + 0.360253^^) 
-I- a^(Afft) (o.00050393 /13 + 2.78092Lq - 2.85654L^ 



0.147051L| 4- 0.865045L^ - 0.165638i| 



0.101931 l: 



Q ' 



(A6) 



where Lq = In ^ and from Eq. (|12p we have = 
8998.080. Our axial-vector hard function for b quarks 



has an extra two-loop singlet piece from the large top- 
bottom mass splitting, H^q^ = H^q + H^q"^^''\ 7J^"gi°t 
was given in Eq. (|13p and involves the real function [sgI ] 

hirt) = 10$(rt)' + 67(rt) + ^ - ^{cisp $(rt)]$(rt) 

+ Cl3[2$(rO]-$(rO'-C(3)}- ;^{2$(rt)Cl2[4$(rt)] 
- 2Cl3[2$(rt)] + Cl3[4$(r,)] + [^jin) + SMnf + C(3)} 



+ \/- - l<j 4(4/i(rO +7(7-t))$(^0 +4Cl2[4$(rt)] 



6$(rt)-6Cl2[2$(rt)] 



Cl2[2$(rt)] +27(rt)$(rt) 



n 



(A7) 



where rt = Q'^/{Am^) and 
$(rt) = arcsin(y7^) , 



7(rt) =ln(2) + -ln(rt), 
Cl2(a;) = Im[Li2(e")] , Cl3(x) - Rc[Li3(e")] , 



ln(2) + -ln(l-r,). 



(A8) 



The resummation of large logs from to is given 
by Uh{Q, tJ'H,IJ'j) in Eq. (jA5|) which is the solution of the 
RGE for the square of the SCET Wilson coefficient [1] 



(A9) 



and the functions u and K are given in Eqs. (jA23p and 
(iA24l) below. 

Finally using results for the convolution of plus- 
functions from Rcf. [ssj we have the momentum space 
formula 



P{Q,k,fij) ^^EP{^,fij,fis) 



E 



n,m,k,l——l 
m+n+l>A; 
A;+l>/ 



Ms 



This result is independent of the dummy variable ^.^^ 
Here £^5^"^ (C, /ij, fJ-s) enco des part o f the running between 
the jet and the soft scale jll2l llisj . 

4'^(C,MJ,A^s)-cxp [2i^(r5 ,js,f^j,fJ's)] (All) 



When convoluted with 5™°'^ we evaluate the right-hand side of 
Eq. ||A10[| for ^ = Qt — 2A{R,fj.s) which simplifies the final 
numerical integration. 
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/ ^\-2iu(rs,tJ.j,tJ.s) exp [2-fEl^iTs,^J'J,^J's)] 



The sum in Eq. ()A10p contains coefficients of tlie mo- 
mentum space soft and jet functions. Sliifting the plus- 
functions so that they have common arguments gives 



' rn=—l 



1 °° C 



(A12) 



n=-l 

Here the thrust soft function coefficients are 



T1=0 



n 



fc=0 



i!fc! 



(A13) 



The soft function is known to 0{al) except for the con- 
stant S3 term [H, [H, US [el 



<S'o(a,s) 



S^i{as) = 1 + 0.349066a^ + (1.26859 -t- 0.0126651 S2)a^ 
+ (1.54284 + 0.00442097 sa + 0.00100786 S3) ^ 
2.07321a^ + (4.80020 - 0.0309077 szja^ , 
-1.69765 as- 6.26659 
- (16.4676 -I- 0.021501 S2)a^ , 
S2{as) = 1.03573 a^- 0.567799 a^, 
Saia,) = 1.44101 a^-f 9.29297 
Siia,) = -1.46525 a^, 

Sr,{as) = -0.611585 al . (A14) 

Note that S2 and S3 are the 0{a^'^) coefficients of the 
non-logarithmic terms in the series expansion of the log- 
arithm of the position space thrust soft function. The 

I 



coefficients appearing in the shifted thrust jet function 
are 



J-i[as,x] =J_i(as) + Jnicts) — — , 

^-^ n + 1 

n=0 

Jn[0is,x] = ^ ^'"LTi,^^' Jn+kias) In'' X , (A15) 



k=0 



and are known up t o 0{a^) except for the constant 
term [20l[58l-l6ll. flTil 



J-i{as) = 1 - 0.608949as - 2.26795a^ 

+ (2.21087 0.00100786 j3) , 

Jo(a^) -0.63662as + 3.00401a^ 
+ 4.45566a^ , 

Ji(a,) = 0.848826as - 0.441765a2 _ u.goSa^ . 

J^las) = -1.0695a2 + 5.36297a^ , 

Jsia,) = 0.360253a2 + 0.169497a^ , 



Ji{a,) = -0.469837a^, 
Jsias) = 0.0764481a^. 

The C distributions are defined as [n > 0] 



(A16) 



9{x) In" X 



da" ^ ^ 



(A17) 



C°Li{x) — £-i{x) = S{x), and for a > —I 



r{x) 



\e^x)^ 




= lim — 




_l_ e-i-o da; 



e{x - e) 



(A18) 



In Eq. (jA10|) we use the coefficients [SS 



d" V{a,b) 
d6" a + b 
n 



6=0 



k db 



kn ■ 



6=0 



a 



1 ' 



fc = -1 , 

< fc < n, 
fc = n + 1 , 



(A19) 
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and the coefficients 



where 



r d™ d" V{a,b) 



da™ d&" a + b 

ni n 
p=0 q=0 



=6=0 



=b=0 



1 1 

+ 



. m + I n + 1 ' 



F(a,6) = 



r(a)r(5) 1 1 



k = -l, 

< k < m + n , 

k = m + n + 1 , 



r{a + b) a b' 

Special cases not covered by the general formulae in Eqs. (|A19[) and (jA20p include 

V-^ (a) = 1 , V,-\a) = a , ¥,7^ (a) = , 



(A20) 



(A21) 



(A22) 



Evolution factors and Anomalous Dimensions 



The evolution factors appearing in Eqs. (jA9|, (|AT0|) . and (|ATT|) are 

/■"^(''^ da 
'^(r, /^o) = 2 / — - r(Q;) 



_ror , as(A^o) /ri /?A l a^(Mo) /ft2 /?2 Ti/? 

3 (4^)3 Lro /?o roV/32 /3oV/32 ^/3o To^J^ M' 



oPo^ 



and 



/7 \ /""^^^^ da 
^2 / A., (mo) p(") 



da^ 



2/32] as(Aio) 



(A23) 



o)V r y Vro ^0^^ ^ 2/3o 47r [Vro/3o /3o'^ 



i?2 In r - 



£2 

Fo ro/3o 



ro/3c 



(1-r) 



(47r)2 



(r^-l) 



rs r2A , S2ri 



^^^^^ ^^^^ B2V.^ln.-^ 



2/3o vro ro/3o 



1 



■ In 7 



(A24) 



where r = as(/i)/as(/io) depends on 4-loop running couplings, and the coefficients are B2 = Pi/I3q — P-i/Pa and 
S3 = -Pl/Pl + 2/3i/32//3g - /33//3o. 



n=0 



n=0 



V47r 



These results are expressed in terms of series expansion 
coefficients of the QCD /3 function /3[as], of r[as] which 
is given by a constant of proportionality times the QCD 
cusp anomalous dimension, and of a non-cusp anomalous The coefficients for Uf = 5 are [6i|,[iil? -[nil 
dimension 7[as], 



n+l 



00 

/3(a.) = -2a.E/3„(^) 



rt+l 



(A25) 



ra=0 



/3o = 23/3 , /3i = 116/3 , ^2 = 180.907, (A26) 
/33 = 4826.16, 
Fg^'P = 16/3, r™"P = 36.8436, T^'^ = 239.208 . 
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For the unknown four-loop cusp anomalous dimension we as{mz) at 4- loops with 5 light flavors we use 
use the Pade approximation assigning 200% uncertainty: 



^ 1 



(A27) 



X 



The anomalous dimensi ons for the ha rd, jet, and soft 
functions are HH, [H [H [ul ITll - fmt 



T^H pCUSp pj ey pC 



ncusp 



Usin) as{mz) 
as{mz) 
X 

X2 



+ 0.401347248 In X (A29) 
[0.01165228 {l-X) + 0.16107961 InX] 



[0.1586117(^2 



0.0599722 {X 



7^ = - 8, 7f = 1.14194, i2 =~ 249.388, 



In X - X"^) + 0.0323244 {(1 - Xf - \n X}] , 



7i 



- 77.3527, -ii 



409.631, 



In 



In In 



where X = 1 + Q!s(m2) ln(^/m2)/3o/(27r) and the dis- 
(A28) played numbers are determined from the Pi in Eq. (jA26p . 

The form in Eq. (jA29[) agrees very well with the numer- 
To determine the strong coupling as{fi) in terms of ical solution of the beta function equation. 

I 



Nonsingular Cross Section Formula 

At 0{al) there is an axial singlet contribution to the nonsingular terms through the threc-parton cut of Fig.[2l which 
is given by the function /singlet appearing in Eq. (|27p for /q"j. The result for this function can be extracted from 
results in Ref. [73 and reads: [n = Q'^ /{Ami)] 



/singlot (''■, 1\ ) 



64 

y 



dyyg{y - l,rt) + (1 - 3r) .g(r, rt) 

l + T ^ 



(A30) 



9{t-, rt) 



2rt 



l — rt T 

rt ■ 



^sin-1 - ^i^sin-i(Vn)] + [sm-\^)]^ ~ [sm-\^)f ~ nlogir) 



4rt(l-r)2 



R-evolution 

Finally we display here the function D^'^) [s^ which ap- 
pears in the solution in Eq. (|4ip of the R-RGE equation 
for A{R, R): 

k 

j=o 

[rhh-j,h)-T{-h-j,t2)], (A31) 

which is real since the complex phase e'^''^ cancels the 
imaginary part coming from the incomplete Gamma 
functions, defined as 



c—l^—x 



(A32) 



Here k is the order of the matrix elements (that is fc = 
for NLL' and NNLL, k = 1 for NNLL' and N^LL, an dfc = 
2 for N^LL'. For lower orders D'^''^ =0). In Eq. (|X3T|) 
we have defined 



72^ 2/?g/?i + Pi - /3o/?2 

•^2 = TTTTTTT TTTT^e. 7l : 



(2/3o)3 



(A33) 



where the R-anomalous dimensions 7/^ were given in 
Eq. (HOD. 

Total Hadronic Cross Section 

The total hadronic QCD cross section, can be evaluated 
in fixed-order perturbation theory with ^ ~ Q, and was 
given in Eq. ([55)1 with the vector QCD results given in 
Eq. (j59p . The function appearing in the singlet contribu- 
tion in Eq. (Ull) at 0{al) is (5^ 



Hn) 



-$(rt)Cl2[2$(rt)]- Cl3[2$(r0]+C3 



1 / 2[l-7(rt)]$(rt)-Cl2[2$(rt)] 

n V n 

-^2Cl2[2<J>(r,)]+2[27(r,)-3]$(r0 



+ 6jirt) + 2^rt)^ ^ 

n 



1 TT^ 



(A34) 



where the necessary functions appear in Eq. (|A8[) . Note 
that we have dropped the four particle cut contribution 
/4 = 7r2/3 — 15/4 since we have not accounted for it in 
the 0{al) nonsingular distribution. 
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Appendix B: Soft Function OPE Matching 

To derive Eq. (|2ip we must demonstrate uniqueness 
of the power correction Q,i and derive its perturbative 
Wilson coefficient to all orders in ag ■ We carry out these 
two parts of the proof in turn. 

Since the operator appearing in the matrix element 
rii is non-local, the proof of uniqueness is more involved 
than for a typical OPE where we could just enumerate 
all local operators of the appropriate dimension. Here we 
are integrating out perturbative soft gluons in Srik,^), 
while retaining nonperturbative soft gluons. The hierar- 
chy between these soft gluons is in their invariant masses, 
fc^ ^ ^QCD- This process can not introduce Wilson lines 
in new light-like directions, nor additional Wilson lines 
following paths in n and n. Thus the Wilson lines will be 
the same as those in the full theory operator, Eq. ([T7|. 
Additional Wilson lines could only be induced by inte- 
grating out collinear or hard gluons, which would yield 
power corrections suppressed by the hard or jet scales. 
The second point to demonstrate is that dimension one 
combinations of derivatives other than id do not lead 
to new nonperturbative matrix elements at this order. 
The key is that for derivative operators inside our vac- 
uum matrix element involving Wilson lines, boost invari- 
ance along the thrust axis relates all matrix elements to 
rii [1^. The proof relics on boost invariance along the 
thrust axis of derivative operators inside the vacuum ma- 
trix element. To see this one defines the transverse energy 
flow operator £t{ii) by its action on states [isl l47j 

Any dimension one derivative operator we might wish to 
consider, such as n ■ d, n ■ d, dt, dz, . . ., or combinations 
thereof, are given by an integral / djy h(r]) £t{v) for an 
appropriate rapidity function h{ri). For example, for the 
thrust derivative id we have h{r]) — e^''''. Boost invari- 
ance implies (65j 

= ^(o|trF^(o)r„(o)*ar„t(o)F;(o)|o) 

= ldv^jO\tTYliO)Yr.iO)£Tiv+ri')YMKiO)\0) 

= ^(0|tr Yl{0)Yr.{0)£T{ri')YMyliO)\0), (B2) 

for arbitrary rj' . The same steps hold for any other deriva- 
tive operator and function and different choices only 
affects the constant calculable prcf actor. This suffices to 
show the second point. 

To derive an all orders expression for the Wilson coeffi- 
cient of r^i we construct an analog of the OPE matching 
done for the soft function in B — >■ Xg^ [1^. The proof 
is considerably simpler for B Xgj because the OPE 
in that case yields local HQET operators. Nevertheless 
the thrust soft function can be manipulated such that a 



similar strategy can be used. Using the thrust axis we 
define hemisphere a where p"*" < p~ and hemisphere b 
where < p'^ . Consider the soft function written as a 
matrix element squared 

Sr{k,fi) = ^^5(fc-fc^-A:,^-)tr (0|F^(0)r„(0)|X) 

^ X 

X (x|y„t(o)F;(o)|o), (B3) 

where the trace is over color, fc"+ = n ■ p'^ is the total 
plus-momentum of the particles in state X in hemisphere 
a and = n ■ p\ is the minus-momenta of particles in 
X in hemisphere b. To carry out the OPE we need to 
consider a state that has overlap with the operator in 
Eq. (|20p . Thus we could replace the vacuum by very soft 
nonperturbative gluons with momenta of ©(Aqcd) and 
then consider matrix elements with perturbative gluons 
having momenta ^ fc ^ Aqcd- Since the OPE is inde- 
pendent of the particular states we choose, we will instead 
consider a simpler alternative in the following. 
First we write the matrix element in Eq. (|B3[) as 

tr(0|F^y„|X)(X|y„tF!|0) (B4) 

= {Q\C,nYlYn^XUnVn){XUnVf,\CnY^YlCn\Q) , 

where C„ and C,n arc non-interacting collinear fields whose 
contractions with the sterile quark u„ and anti-quark Vn 
are chosen with a normalization to reproduce the origi- 
nal matrix element (and a sum over their color correctly 
reproduces the trace). Here it„ should be thought of as a 
very energetic collinear quark in hemisphere a with large 
label momentum p~ , and zero residual momentum. The 
large momentum is conserved by soft interactions from 
the Wilson lines due to the SCET multipole expansion. 
Here the plus-momentum of m„ is included into , but 
is zero and does not contribute to the (5-function. The 
same is true for Vn which has zero minus-momentum, 
large label momentum, and is always in hemisphere 
b. We introduced w„ and Vn so that we can use them 
to systematically add a very soft momentum to the end 
of the Wilson lines (at oo). They provide a convenient 
state with which to carry out the OPE, because there is 
nonzero overlap taking only the 1 out of the Wilson lines, 
Y . In particular they allow us to perform the OPE and 
pick out the id present in f^i at tree level, without the 
necessity to add explicit soft gluons with momenta ^ k. 

To carry out the OPE we now give m„ a very small 
soft momentum ^+ and Vn a very small soft momentum 
and denote them by ufj and v^^ respectively. These 
particles are kept on-shell by adjusting their large label 
_L-momenta so that = Pn±/Pn and = p'^±/pt- 
Due to the multipole expansion these _L-momcnta have 
no influence on diagrams with perturbative soft gluons 
having momenta k <^ pn±_ — —pn±- The Wilson line 
propagators reduce to the same as before, such as 



k+pn + i+Pn + pIj^ k+ ■ 
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FIG. 21: Amplitudes for zero and one soft gluon. 



This property is familiar in SCET where soft couphngs 
to energetic colhnear quarks in SCET remain eikonal for 
any values of the quark's large momenta by using the 
equations of motion, as long as the final particles are on- 
shell. Thus at any order in perturbation theory, with any 
number of soft gluons and soft quarks of momenta ~ k 
in the matrix elements, the only change caused by £^ is 
on the S{k - fc°+ - fc^") in Eq. ((B3t which is shifted to 
6{k - i - - fcj-), where £ = £+ + £-. Expanding 
with £ <^ k the matrix element with this choice of state 
evaluates to 

{k-£,fi) = 5P"* (fc) - '^'^"■7^^^ £+.... (B6) 

At lowest order in emission of very soft gluons ~ Aqcd 
the corresponding matrix element in the lower energy 
theory is 

^{OKnYlY„Cn^^\uivi){uivi\CnY^YlCf^\0) ^£. 

(B7) 

Virtual radiative corrections do not correct this result 
since they are scaleless and vanish in pure dimensional 
regularization. Thus we can identify £ — >■ 2rii in 
Eq. (jB6|) . and this then yields the stated result for the 
OPE in Eq. ((2T|) . 



Appendix C: Operator Expansion for the First 
Thrust Moment 



For moment integrals of the thrust distribution over 
r e [0,1/2] there is not a hierarchy of scales that in- 
duces large logs, and one may formulate the theoreti- 
cal result in terms of an expansion in and Aqcd/Q- 
The zero'th moment of thrust is just the total cross sec- 
tion for e+e^ — > hadrons, and the power corr ectio ns are 
formulated in terms of the well known OPE jl22l |. For 
higher moments the fact that thrust constrains a non- 
trivial combination of final state momenta makes carry- 
ing out an OPE more difficult. For example, when we 



weigh the integral by a power of thrust it is not possible 
to collapse all propagators to a point, so the nonpertur- 
bative parameters are no longer constrained to be given 
by a basis of local operators. In the effective coupling 
model [ssf the same nonperturbativc parameter ao that 
appears for the thrust distribution, also occurs in the first 
moment. However it is not clear to what level of accuracy 
this carries over to a field theoretical description of power 
corrections derived from QCD. In this appendix we show 
how one can carry out an OPE for the 1st moment of 
the thrust distribution, and demonstrate that at leading 
order it only involves the same nonperturbativc matrix 
element fti from Eq. 

To carry out an OPE for the thrust moment we can 
work order by order in the hard as{Q) expansion, and an- 
alyze direct computations where we couple soft nonper- 
turbativc gluons to hard partons in Feynman diagrams. 
The appropriate non-local operator(s) appearing in the 
expansion will be identified by the structure of the am- 
plitudes in this computation. In the following discussion 
the soft gluons will not be treated as final state particles 
for which there is a phase space integral, but rather as 
a means of probing the structure of the nonperturbativc 
operator. The lowest order graphs with zero or one soft 
gluon and a virtual photon current (for simplicity) are 
shown in Fig. [5TJ Here fc'^ ^ Aqcd is soft, and p'' ^ Q, 
p''^ ^ Q are hard momenta. To carry out the OPE we cal- 
culate and square the on-shell amplitude, A^^AIJ^, where 
/i, V are the virtual photon current indices. We sum over 
the final quark/antiquark spins since these particles are 
hard and are being integrated out. On the other hand 
the gluon vector indices a, a' are left uncontracted and 
are used to help in identifying the operator for the non- 
perturbativc matrix element. For simplicity, the indices 
a and a' are suppressed in writing down the amplitudes 
below. We start out without making restrictions on the 
number of gluons coming from and A^^*, which cor- 
responds to directly matching onto the nonperturbativc 
operator, without considering the final vacuum matrix 
element which gives a nonperturbativc parameter. Since 
p2 = p'2 _ ^2 _ Q ^j-^g denominators of the propagators 
in the one gluon graphs reduce to 2p ■ k and 2p' ■ k. In 
the numerators we can drop |^'s relative to the large ^ 
and The interference between the zero and one gluon 
amplitudes gives 



2gT 



A r 



Nr. 



p" 



p' ■ k p ■ k 



(CI) 



The interference with one gluon from each of and 



Mr: is 



7Vctr[|i7^/7^1 



(C2) 



Nr. 



plap,a 



{p"p" 



■ p'^'p"' 



{p-kY {p' -ky {p-k)(p' -k) 



Continuing in this fashion with any number of gluons 
from ^-^^id any number from Ai^C we always find 
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the tree level amplitude squared with no soft gluons, 
Nctrl^^i^i)'^"], times an amplitude from the soft gluons. 

Since the hard quarks are on-shell and back-to-back 
their four-momenta are given by light-like vectors along 
the thrust axis, 



P 



ifj, 



,n-p 



(C3) 



in J dT{T / a){da / dr) is given by 



. \2p-k 2p'-k 
T ~ mm 



1 . 

—mm 

Q 



n ■ k,n ■ k 



(C4) 



— |n • kO{n ■ k — n ■ k) + n ■ kO{n • A; — n • fc)| , 



up to power corrections beyond those considered here. 
Here = (l,t) and = (1,— t) are identical to the 
n and n appearing in Eq. (|6]). Using Eq. (|C3|) the soft 
gluon amplitudes in Eqs. (jCip and (jC2l) are cikonal with 

precisely the right factors to come from the Y^^ (0), Yn{0), 
Fj(0), F^(0) in the f^i matrix element in Eq. ([6|). 

For the first moment observable we can focus on am- 
plitudes that have the same number of gluons in A^^ 

and M.'{*, and at least one gluon for the id operation in 
Eq. ([6]) to act on. Since the gluon is soft, the factor of r 



and is exactly equal to id given in Eq. ([7]) acting on 
the soft gluon in Fig. [5TJ Hence in the first moment 
of thrust we find that r together with the soft gluon 
amplitude give precisely 2Q,i/Q^ with the vacuum matrix 
in Eq. ([6]) (where the trace comes from the sum over color 
for the final state quarks). The remaining iVc tr[/57^/5'7''] 
amplitude goes together with the two-body phase space 
to yield the tree level cross section ctq. Together these 
results yield Eq. p5|) for the lowest order OPE for the 
first moment of thrust. 
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